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The  primary  goal  of  our  work  was  to  provide  accurate  information  about  velocity  and 
density  field  in  the  near  wake  region  of  a  flow  around  a  sphere  in  a  stratified  fluid.  This  is  a 
canonical  configuration  used  in  investigating  a  structure  of  stratified  wakes  in  experiments 
(e.g.,  Spedding  et  al.  1997,  2001,  2002)  and  in  numerical  simulations  (c.g.,  Dommermuth  ct 
al.  (2002),  Diainessis  et  al.  (2005),  Brueker  and  Sarkar  (2009)).  One  of  the  difficulties  in 
comparing  results  of  numerical  simulations  with  experiments  is  lack  of  detailed  information 
about  near  wake  that  is  needed  to  initialize  simulations.  Despite  importance  of  this  config¬ 
uration  it  appears  that  only  one  numerical  investigation  has  been  performed  that  contains 
an  actual  sphere  in  the  computational  domain  (Hanazaki,  1988).  However,  the  simulated 
flow  in  that  work  is  at  low  Reynolds  number  (Re=200),  well  below  turbulent  regime.  Our 
research  work  was  designed  to  fill  this  gap  and  to  provide  additional  information  useful  for 
numerical  simulations,  turbulence  modelling,  and  the  interpretation  of  some  experimental 
techniques. 


1  Summary  of  Research  Activities 


The  computations  in  the  project  were  supervised  by  Prof.  J.A.  Domaradzki.  A  PhD  student 
in  the  AME  department,  Trevor  Orr,  has  worked  on  the  project  full  time;  this  is  the  topic 
of  his  Ph.D.  Thesis,  with  Qualifying  Exam  passed  in  September  2010.  We  have  established 
collaboration  with  Prof.  George  Constantinescu  and  acquired  a  serial  code  for  simulating 
a  flow  around  a  sphere  in  curvilinear  coordinates  using  detached  eddy  simulations  (DES) 
approach.  The  code  was  first  implemented  and  tested  for  the  flow  without  stratification. 
Subsequently,  it  has  been  modified  to  include  effects  of  stratification,  again  in  the  framework 
of  DES.  Since  the  modified  code  already  contained  turbulence  model  we  were  able  to  move 
directly  to  simulating  flows  at  experimental  Reynolds  numbers  where  comparisons  with  a 
single  existing  dataset  of  Spedding  for  a  near  wake  are  possible  (Re=5000).  The  numerical 
results  for  the  mean  wake  velocity  were  found  to  be  in  good  agreement  with  experiments  but 
rms  velocities  showed  poor  agreement.  The  initial  hypothesis  was  that  observed  results  for 
rms  velocities  were  caused  by  internal  waves  reflecting  from  the  computational  boundaries. 
To  decrease  this  effect  various  numerical  absorbing  layers  were  implemented  and  tested  but  no 
significant  improvement  has  been  observed.  We  undertook  a  systematic  analysis  of  possible 
origins  of  the  rms  turbulent  kinetic  energy  build  up  in  the  DES  simulations  for  stratified 
wake:  is  it  caused  by  numerical  resolution?;  by  the  turbulence  model?;  by  the  boundary 
conditions?;  by  stratification  and  internal  waves?  To  increase  the  numerical  resolution  a 
parallel  version  of  the  code  has  been  implemented  on  the  USC  Linux  cluster.  To  determine 
the  influence  of  the  model  we  have  initiated  DNS  runs  (without  a  model)  for  a  stratified 
case.  To  determine  the  effects  of  boundary  conditions  we  have  initiated  DNS  and  DES  for 
an  unstratified  case.  The  results  for  stratified  flow  cases  for  Re— 200  and  Fd—  0.125  and  0.250 
compare  very  well  with  Hanazaki  (1988);  for  unstratified  flow  cases  for  Re=300  we  obtain  a 
good  comparison  with  Johnson  and  Patel  (1999).  All  relevant  results  are  summarized  in  a 
separate  section  below. 

During  the  duration  of  the  project  PI  was  involved  in  additional,  related  research  with 
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unsupported  Ph.D.  student  on  subject  of  LES  modeling.  The  relevant  paper  (attached) 
acknowledges  partial  ONR  support  for  that  work.  Another  invited  review  paper,  currently 
in  print,  on  subject  of  LES  modeling  without  eddy  viscosity  models  also  acknowledges  ONR 
support  and  is  attached. 
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4  Research  Results 
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tailed  information  about  research  results  summarized  above. 
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Chapter  1 

Introduction  and  Background 


1.1  Introduction 

1.1.1  A  Description  of  the  Physical  Problem 

The  physical  problem  involves  the  numerical  modeling  of  a  sphere  traveling  horizontally  through 
a  stable,  linearly  stratified  fluid  of  constant  density  gradient  at  moderate  Reynolds  (Re)  numbers 
for  a  range  of  Fronde  (Fd)  numbers.  The  investigation  will  study  the  effects  seen  by  the  presence  of 
stratification  in  both  laminar  (low  to  moderate  Re)  and  turbulent  (moderate  to  high  Re)  regimes, 
but  sphere  motion  at  very  low  speeds  (creeping  flow)  is  not  included  in  this  work. 

Investigation  of  behavior  of  the  flow  around  a  body  in  a  stratified  fluid  is  of  interest  to  engi¬ 
neering  applications  as  well  as  some  flows  at  the  geophysical  scale;  flows  over  mountain  ridges  for 
instance.  Engineering  applications  may  involve  flow  around  a  submarine.  Because  the  applications 
of  these  flows  are  quite  varied,  there  are  so-called  canonical  geometries  which  aid  computational 
modelers  in  verifying  the  behavior  of  the  computational  models.  The  sphere  is  generally  accepted 
to  be  one  of  these  canonical  problems  because  of  its  simple  geometric  shape  as  it  contains  a  perfect 
symmetry  about  its  center  point. 

Flow  around  a  submersed  object  in  a  stably  stratified  fluid  has  several  physical  features  that  are 
not  found  in  a  similar  flow  without  density  stratification;  i  e.  a  homogeneous  flow.  This  is  because 
of  restorative  buoyancy  forces  within  the  fluid  that  return  any  disturbances  in  the  stratification 
profile  to  its  previously  undisturbed,  "at  rest,”  condition. 


1.1.2  The  Governing  Equations 
Equations  of  motion  and  assumptions 

The  problem  of  the  sphere  traveling  horizontally  through  a  continuously  stratified  fluid  is  gov¬ 
erned  by  a  particular  set  of  the  Navier-Stokes  equations  of  motion  and  the  equation  of  continuity: 


dpux  dpu%Uj  OP  d2ux 

dt  +  dxj  dxx  ^ Ox j2 


0,  fi  =  P9$n 


dp  dpuj 
dt  dxx 


(1.1) 

(1.2) 
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where  e*  represents  the  orthnormal  coordinate  vectors  in  a  Cartesian  coordinate  system,  u  =  u%et 
is  the  velocity  vector  field,  P  is  pressure,  p  is  the  density,  and  ft  is  the  component  of  body  forcing 
due  to  gravity  acting  on  individual  fluid  particles  per  unit  volume. 

It  is  desirable  to  decompose  both  the  pressure  and  density  field: 

p  (x ,  t)  =  po  +  p  (z)  +  p  (x,  t)  (1.3) 

P(f,  t)  =  P0  +  P(z)  +p  {x,t)  (1.4) 

This  decomposition  uses  the  concept  of  hydrostatic  balance  between  a  background  pressure,  p{z), 
and  a  background  density  distribution,  p(z).  Physically,  this  translates  into  a  balance  between  the 
gravity  force  acting  on  a  fluid  particle  and  the  restoring  buoyancy  forces  exerted  by  any  surrounding 
fluid  particles  while  the  fluid  is  at  rest  such  that  only  disturbances  to  this  equilibrium  exhibit  any 
net  effect.  Mathematically,  the  hydrostatic  balance  aids  in  the  proper  construction  of  boundary 
conditions,  especially  in  fluids  of  infinite  depth.  Without  the  hydrostatic  balance,  the  pressure 
conditions  must  be  explicitly  set  at  the  boundaries.  The  primed  quantities  P  and  p  are  referred 
to  as  the  perturbation  pressure  and  perturbation  density,  respectively. 

The  equation  of  continuity,  Eq.  1.2,  enforces  the  conservation  of  mass.  The  Boussinessq  ap¬ 
proximation  is  used  [26].  This  assumes  that  any  changes  in  density  are  only  accounted  for  in  the 
body  forcing  term.  In  addition,  any  effect  of  a  fluid  particles'  density  change  on  the  equation 
of  continuity  is  negligible  and  leaves  the  continuity  equation  as  stating  that  the  velocity  field  is 
divergence  free. 

Substituting  the  two  ideas  above  into  the  Eqs.  1.1  and  1.2: 


dut  diii  dp  p  d2ut 

dt  Uj  dxj  ^  dxi  p^dxj2 


(1.5) 


dui 

dxi 


=  0 


(1.6) 


An  equation  must  also  be  presented  for  evolution  of  the  perturbation  density  quantity.  As 
stated  by  Kundu,  an  equation  for  the  description  of  the  density  perturbations  is  provided  by  the 
concept  of  incompressibility  (e.g.  Kundu  pg.  248  [26]);  with  the  current  formulation  containing 
dissipative  effects.  The  equation  for  the  convection  of  density  perturbations,  known  as  "the  density 
equation”  becomes: 


dp_ 

dt 


+  u3-k 


o2p 

dxj 2 


=  0 


(1.7) 


where  k  is  the  molecular  diffusion  coefficient. 


Non-dimensionalization  of  the  governing  equations 

It  is  beneficial  to  complete  a  non-dimensionalization  of  the  governing  equations  to  distill  the 
physics  of  the  problem  into  a  problem  controlled  by  parameters  of  similtude.  Non-dimensionalized 
quantities  within  this  section  will  be  denoted  by  an  asterik  (*);  which  will  be  dropped  after  this 
section  for  aesthetics  of  the  recast  equations.  The  diameter  of  the  sphere,  D,  the  speed  of  the 
sphere,  U,  the  characteristic  density,  po,  and  the  background  density  stratification,  dp(z)/dz ,  are 
intuitive  choices  for  the  characteristic  quantities  of  the  problem.  Recasting  Eqs.  (1.5)  and  (1.6)  in 
non-dimensional  form: 
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e3  =  0 


(1.8) 


dux*  +  dut*  dp 
dt*  J  dxj*  +  dxx* 


1  d2ux* 
Re  dxj  * 2 


/  * 

P 

F<P 


dut * 
dxx * 


=  0 


d£_  ,dp^_  .  i 

dt *  dxj*  U3  Pedxj *2 


=  0 


(1.9) 

(1.10) 


where  f*  =  [/*/£>,  x*  =  x/D ,  iT  =  u/C/,  p*  =  p'/poU2,  p*  =  p  /  [-D  dp{z)/dz] ,  Re  - 
DU jv,  Fd  =  U/(N  D),  5c  =  «/(/,  and  Pe  =  1/  5ci?e  have  been  used  as  non-dimensionalization 
parameters. 

Here,  iV(z)2  =  [— dp(z)/dz]  <///>o  is  the  Brunt-Vaisala  frequency.  In  this  study,  the  density 
stratification  is  continuous  and  linear  in  undisturbed  conditions,  resulting  in  N(z)  =  N  =  constant, 
but  no  generality  has  been  lost  in  the  governing  equations.  Through  non-dimensionalization,  if  the 
Schmidt  number  is  taken  to  be  unity,  it  is  apparent  that  Eqs.  (1.8)-(1.10)  are  controlled  by  the 
parameter  of  the  Froude  and  Reynolds  number. 


1.1.3  Stratified  Fluids 

The  fluid  may  be  a  liquid  or  gas  of  a  varying  density  where  the  continuum  approximation  is  a 
valid  assumption.  If  an  infinite-in-tho-horizontal  volume  of  stationary  fluid  is  positioned  directly 
above  a  second  similarly  constructed  volume  of  fluid,  and  the  density  of  the  fluid  in  the  first  volume 
is  less  than  that  contained  in  the  second,  and  this  pattern  is  repeated  for  all  layers  in  the  vertical, 
then  the  fluid  may  be  considered  stably  stratified.  A  stable,  continuously  stratified  fluid  is  the 
occurrence  where  these  stacked  layers  become  infinitesimally  thill  in  the  vertical  direction,  and  the 
density  gradient  of  the  stratification  maintains  the  stable  condition  for  all  these  layers. 

The  density  gradient  of  the  stratification  is  embodied  within  the  so-called  Brunt-Vaisala  fre¬ 
quency  which  is  commonly  denoted  as  V(z).  Physically,  the  Brunt-Vaisala  frequency  represents 
the  frequency  at  which  a  fluid  particle  that  has  been  vertically  displaced  by  some  amount  will 
oscillate  about  its  original  equilibrium  state.  As  the  stratification  gradient  of  the  continuously 
stratified  fluid  becomes  increasingly  severe,  the  Brunt-Vaisala  frequency  increases  proportionally 
to  the  value  of  the  local  density  gradient.  For  the  steady  motions  of  submersed  bodies,  the  Froude 
number  becomes  the  analogue  of  the  Brunt-Vaisala  frequency  ,  and  as  the  Froude  number  decreases, 
the  stratification  of  the  fluid  becomes  increasingly  severe. 

In  low  enough  Froude  number  flows,  a  unique  effect  of  stratification  (excluding  creeping  flow)  is 
referred  to  as  the  11  blocking  phenomenon.”  The  blocking  phenomenon  is  not  a  primary  focus  of  this 
work  but  may  be  encountered  for  some  flows  contained  within  this  investigation  and  mention  of  the 
physics  of  these  flows  is  prudent.  The  blocking  phenomenon  is  a  result  of  the  flow  lacking  enough 
kinetic  energy  to  convert  into  the  potential  energy  required  to  move  vertically  over  a  submersed 
obstacle.  In  two-dimensional  flows,  the  effect  of  this  energy  deficit  is  that  fluid  particles  are  trapped 
on  the  upwind  side  of  the  object  and  are  forced  to  move  in  unison  with  the  obstacle’s  motion.  The 
fluid  behaving  in  this  fashion  is  referred  to  as  ’’blocked.”  In  three-dimensional  flows,  the  fluid  has 
the  ability  to  flow  around  the  sides  of  the  object  if  there  is  not  enough  kinetic  energy  to  overcome 
the  potential  energy  required  to  move  over  the  object,  but  even  in  this  case  some  blocking  effects 
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may  still  be  seen  in  some  planes  of  symmetry.  A  review  article  by  Baines  in  1987  [2]  provides  a 
review  into  the  topic  of  the  blocking  phenomena. 

As  the  Froude  number  increases,  the  effects  of  stratification  become  less  apparent,  the  Brunt- 
Vaisala  frequency  decreases,  and  as  the  Froude  number  approaches  the  infinite  limit,  the  flow 
behaves  as  a  homogeneous  fluid  flow. 

1.2  Stratified  Flow  Around  a  Sphere 

A  submersed  body  is  assumed  to  be  in  a  steady  lion-oscillatory  motion.  References  to  the  flow 
around  the  body  indicate  that  the  observer’s  reference  frame  is  fixed  with  respect  to  the  body  as 
it  moves  through  a  stationary  fluid,  and  variances  to  this  convention  are  otherwise  noted.  Ref¬ 
erences  to  ”the  sphere”  are  intended  to  reference  the  sphere  traveling  horizontally  through  the 
linearly  stratified  fluid  at  a  constant  speed.  Flows  that  maintain  an  orderly,  predictable  behavior 
insensitive  to  perturbations  or  initial  conditions  are  called  laminar.  Flows  that  do  not  maintain 
these  properties  are  generally  referred  to  as  turbulent;  although  there  is  no  universally  accepted 
definition  of  turbulence.  A  sketch  of  the  physical  problem  is  seen  in  Figure  4.1. 

1.2.1  Behavior  of  Flow  Around  the  Sphere 

The  first  region  of  physical  interest  is  the  behavior  of  the  flow  around  the  body  including  the 
flow  close  to  the  surface  of  the  sphere.  Early  focus  on  the  flow  upwind  and  around  a  submersed 
obstacle  in  a  stratified  fluid  is  typified  by  the  flow  over  mountains  or  hill-shaped  objects,  especially 
in  some  of  the  earliest  work  on  the  topic.  Although  physically  focused  on  a  geophysical  problem, 
these  studies  are  typically  intended  to  describe  the  behavior  of  fluid  particles  as  they  approach 
and  eventually  attempt  to  flow  around  a  submersed  obstacle,  which  is  analogous  to  the  flow  of 
fluid  particles  around  the  sphere.  Attempts  to  quantify  the  ability  of  fluid  particles  to  traverse 
an  obstacle  is  discussed  by  Sheppard  1956  [43  by  considering  energy  arguments.  The  efforts  of 
Drazin  1961  [17]  attempt  to  resolve  the  question  of  whether  a  fluid  particle  will  flow  over  or  around 
a  three-dimensional  object.  A  review  of  the  various  attempts  to  quantify  and  refine  the  theories 
presented  by  [43]  and  [17]  are  discussed  in  Hanazaki  1988  [21].  The  relevant  summaries  to  these 
works  as  applied  to  the  sphere  case  is  that  fluid  particles  in  the  region  affected  by  the  presence 
of  the  sphere  gain  a  directional  preference  as  to  whether  they  will  flow  over  or  around  the  side  of 
the  sphere  because  of  the  addition  of  the  buoyancy  forcing  due  to  gravity.  The  further  away  from 
the  sphere’s  region  of  influence,  the  less  this  preferred-direction  effect  is  seen  and  the  free  stream 
motions  of  the  fluid  are  generally  considered  to  be  unaffected.  Thus,  the  flow  around  the  sphere  is 
mainly  concerned  with  particle  paths.  The  behavior  of  these  particle  paths  becomes  more  effected 
by  the  presence  of  the  sphere  the  closer  the  flow  advects  the  particles  to  the  sphere  surface. 

Within  this  region  close  to  the  sphere,  the  fluid  particles  on  the  surface  of  the  sphere  are  assumed 
to  be  ’’stuck”  to  the  surface  and  cannot  be  advected  away  from  the  surface.  This  is  referred  to  as  the 
”  no-slip  condition,”  and  is  used  for  both  homogeneous  and  stratified  flows.  The  no-slip  condition 
forces  the  relative  velocity  of  the  attached  fluid  to  be  zero  with  respect  to  the  surrounding  free 
stream  flow.  Viscous  forces  develop  within  the  fluid  due  to  this  momentum  mismatch  near  the 
surface  of  the  sphere.  This  region  is  referred  to  as  the  boundary  layer  [40].  The  Reynolds  number 
is  the  ratio  of  momentum  to  viscous  effects  and  serves  as  an  indicator  of  the  behavior  of  fluid  flow 
near  physical  boundaries  and  subsequent  flow  development  downstream  of  these  boundaries. 
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The  boundary  layer  of  the  sphere  exhibits  different  characteristics  depending  upon  the  Reynolds 
number,  the  Froude  number,  and  possibly  the  diffusivity  [22]  of  the  flow.  In  turn,  the  engineering 
characteristics  (e.g.  drag  or  shedding  cycles)  of  the  sphere  may  depend  on  the  characteristics  of 
the  boundary  layer.  The  boundary  layer  on  the  surface  of  the  sphere  may  be  roughly  characterized 
by  two  regimes  separated  by  a  ’’critical”  Reynolds  (Re^t)  number  [26].  In  this  study,  focus  will 
be  maintained  in  the  ’’sub-critical”  regime,  Re  <  where  the  boundary  layer  flow  is  laminar. 

In  a  homogeneous  fluid,  viscosity  affects  the  flow  around  the  sphere  in  the  laminar  boundary 
flow  around  the  sphere  below  a  Reynolds  number  of  approximately  25  such  that  the  particles 
close  to  the  surface  of  the  sphere  travel  a  path  parallel  to  the  contour  of  the  sphere  surface  in 
an  axisymmetric  manner  [52].  Above  a  Reynolds  number  of  25,  and  below  a  Reynolds  number 
of  roughly  200  [25],  the  flow  behaves  in  a  steady,  axisymmetric  manner,  but  the  flow  creates  a 
so-called  recirculation  bubble,  or  vortex  ring,  that  forms  behind  the  sphere  [52].  In  this  case,  the 
convecting  fluid  particles  within  the  flow  do  not  follow  the  contour  of  the  sphere’s  surface.  Instead, 
the  particles  are  convected  away  from  the  surface  of  the  sphere,  and  the  flow  velocity  field  creates 
a  stationary  region  of  vorticity  behind  the  sphere  known  as  a  ’’separation  bubble.”  The  point  at 
which  the  fluid  particles  traveling  within  the  boundary  layer  stop  following  the  surface  contour  is 
commonly  referred  to  as  the  ’’separation  point,”  and  the  boundary  layer  after  the  separation  point 
itself  is  often  referred  to  as  being  ’’separated.” 

These  two  particular  cases  of  boundary  layer  behavior  dependent  on  the  Reynolds  number 
provide  an  opportunity  to  discuss  separation  in  the  context  of  boundary  layer  flows.  Qualitatively, 
separation  is  the  terminology  given  to  the  process  of  the  fluid  particle  within  the  boundary  layer 
erupting  into  the  free  stream  [54].  The  exact  definition  of  a  separation  point  is  continually  evolving 
within  literature  starting  with  the  original  definition  of  Prandtl  [39]  with  additional  modifications 
to  account  for  behavior  of  unsteady  flows  (e.g.  [42]  or  [20]).  These  modified  definitions  all  tend 
to  collapse  back  to  the  original  Prandtl  definition  for  the  case  of  steady,  or  near-steady  flows. 
Thus,  in  this  study,  the  choice  for  the  definition  of  the  ’’separation  point”  is  where  the  surface  skin 
friction  coefficient  initially  approaches  zero,  and  a  ’’separation  line”  is  the  interpolation  between  the 
individual  separation  points.  Boundary  layers  are  then  considered  ’’separated”  after  the  separation 
point,  except  in  the  rare  exception  of  reversed  boundary  layer  flows  that  do  not  qualitatively 
coincide  with  the  physical  description  of  [54].  The  separation  behavior  of  the  sphere  traveling 
through  the  stratified  fluid  is  of  interest,  especially  as  it  relates  to  the  separation  behavior  seen  in 
a  homogeneous  fluid. 

Claims  to  the  first  study  of  the  relationship  between  the  separation  line  on  the  sphere  in  con¬ 
tinuously  stratified  fluid  arrives  from  Lofquist  &  Purtell  1984  [34].  There  is  additional  work  done 
by  Sysoeva  &  Chashechkin  1986  [51];  and  qualitative  results  again  by  Chashechkin  &  Sysoeva  in 
1988  [8],  specifically  at  that  time  with  Fd  <  0.2.  The  first  computational  work  done  in  the  area  of 
the  stratified  fluid’s  effect  on  separation  is  from  Hanazaki  1988  [21].  The  result  of  these  investiga¬ 
tions  is  that  separation  behavior  in  the  vertical  plane  is  delayed  as  compared  to  the  homogeneous 
case  as  the  separation  point  moves  further  and  further  downstream  on  the  surface  of  the  sphere. 
However  as  the  severity  of  the  stratification  increases,  the  vertical  separation  point  tends  to  move 
slightly  upwind  as  the  flow  approaches  a  quasi  two-dimensional  regime  induced  by  the  extreme 
stratification. 

According  to  Hanazaki  1988,  at  Re  =  200,  Fd  =  0.5,  separation  in  both  vertical  and  horizontal 
planes  is  totally  suppressed  although  Sysoeva  &  Chashechkin  1986  do  not  see  a  total  suppression 
of  separation  at  that  Reynolds  number.  Lofquist  h  Purtell  1984  do  not  register  a  total  suppression 
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of  separation  at  700  <  Re  <  1000,  Fd  =  0.5  either  but  do  measure  a  maximum  separation  angle 
of  160  degrees.  Additional  work  by  Chomaz  et  al.  in  1992  [10]  present  separation  studies  at  both 
low  and  high  Reynolds  numbers.  Their  experimental  results  indicate  that  a  difference  in  Reynolds 
number  accounts  for  the  discrepancy  of  separation  suppression  between  [34]  and  [21].  They  also 
discount  the  discrepancy  at  Re  =  200  between  Hanazaki  and  Sysoeva  &  Chashechkin  as  likely  due 
to  the  apparent  inability  of  the  shadowgraph  technique  to  accurately  capture  separation  location 
[10]. 

The  sketch  of  Figure  5  in  Chomaz  et  al  [10]  extends  the  earlier  two  dimensional  outlook  from 
the  vertical  and  horizontal  planes  and  is  a  qualitative  interpretation  of  how  stratification  affects 
the  separation  line  on  the  lee  side  of  the  sphere.  Generally,  as  the  Reynolds  number  is  decreased, 
the  surface  area  in  the  separated  region  decreases  [10].  There  have  been  no  computational  studies 
of  the  effects  of  stratification  on  the  shape  of  the  separation  line  for  the  sphere  at  any  Reynolds 
numbers  other  than  what  may  be  implied  from  the  data  from  Hanazaki  1988  at  a  Revnolds  number 
of  200. 

1.2.2  The  Wake  Region 

The  second  physical  region  of  interest  is  the  wake  region  located  on  the  downwind,  lee,  side  of  the 
sphere.  There  are  several  types  of  laminar  wake  behaviors  behind  the  sphere;  depending  generally 
on  the  value  of  the  Reynolds  number  and  Froude  number.  The  various  types  of  wakes  found 
in  homogeneous  fluids  will  be  listed  as:  steady  asymmetric  vortex,  unsteady  symmetric  vortex 
shedding,  asymmetric  vortex  shedding,  but  reviews  of  the  Reynolds  numbers  these  occur  at  may  be 
found  in  [25]  or  [10].  The  properties  of  the  wake  may  change  depending  on  Froude  number  as  well, 
but  most  structures  present  in  the  stratified  wake  are  identifiable  as  some  form  of  the  previously 
mentioned  homogeneous  wake  behaviors.  One  unique  regime  found  in  the  stratified  fluid  flows  is 
the  existence  of  the  standing  lee  wave,  which  is  a  laminar  mechanism.  Above  a  Reynolds  number  of 
2000,  it  is  expected  that  the  wake  region  will  become  turbulent  for  all  Froude  numbers  higher  than 
2.  Below  a  Froude  number  of  2,  the  balance  between  strength  of  stratification  and  inertial  effects 
(suggested  as  Fd  >  C  *  Re~ 1//2  by  [31]  where  C  is  a  calibrated  constant)  will  determine  whether 
the  wake  is  turbulent  or  laminar  according  to  Figures  2  &  3  in  Lin  et  al.  [32];  reproduced  here  for 
convenience  as  Figure  4.5. 

Theoretical  work  done  in  the  wake  of  the  sphere  utilizes  linear  theory  in  the  context  of  a  moving 
point  source  that  is  naturally  generating  internal  waves  (e.g.  Lighthill  [28]).  Internal  waves  may 
exist  for  all  stratified  flows.  For  flows  involving  a  density  interface  of  layered  fluids,  the  internal 
waves  are  concentrated  along  the  density  interface  between  layers  whereas  for  continuously  stratified 
flows,  internal  waves  may  propagate  freely  throughout  the  fluid.  The  result  of  linear  theory  gives 
insight  into  the  behavior  of  the  internal  waves.  Linear  analyses  are  generally  presented  in  terms  of 
the  lines  of  constant  phase  of  any  waves  that  propagate  outward  from  the  sphere.  Linear  theory 
may  also  be  used  as  a  reference  to  measure  the  wavelength  of  any  waves  present  behind  the  sphere. 
Because  of  this,  comparisons  of  experimental  and  computational  observations  are  often  made  in 
some  fashion  with  linear  theory. 

As  mentioned,  a  unique  characteristic  of  the  flow  in  a  stratified  fluid  is  the  existence  of  dominant 
lee  waves.  A  stationary  lee  wave  is  an  internal  wave  propagating  with  a  phase  speed  relative  to 
that  of  the  sphere  it  follows  such  that  in  the  frame  of  reference  to  the  body,  it  is  stationary.  The 
existence  of  a  standing  lee  wave  is  seen  in  the  work  of  Hanazaki  1988.  Hanazaki  1988  [21]  conducted 
simulations  at  a  steady  Reynolds  number  of  200  for  a  range  of  Froude  number  of  0.15  to  200,  which 
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is  effectively  an  unstratified  flow  at  that  Reynolds  number.  Hanazaki’s  work  makes  comprehensive 
numerical  study  of  the  relationship  between  the  drag  on  the  sphere  by  standing  lee  wave  behavior 
in  the  near  wake  as  well  as  provides  qualitative  information  about  the  flow  field  surrounding  a 
sphere  in  a  continuously  stratified  fluid.  The  conclusion  of  the  Hanazaki  work  is  that  the  linear 
theory  predictions  are  generally  well  represented  [21]. 

The  work  of  Lofquist  &  Purtell  in  1984  [34]  also  investigates  the  drag  coefficient  of  a  sphere 
moving  through  a  linearly  stratified  fluid.  Lofquist  &  Purtell  1984  suggest  that  regardless  of 
Reynolds  number,  the  lee  wave,  which  is  dependent  on  the  stratification,  induces  drag  [34].  It  is 
suggested  that  there  is  an  increase  in  drag  coefficient  that  is  independent  of  Reynolds  number  as 
the  fluid  becomes  increasingly  stratified  until  a  Froude  number  of  about  0.125.  Below  that  Froude 
number,  the  drag  coefficient  begins  to  decrease  as  Froude  number  decreases.  This  is  shown  in 
Figure  4  of  that  work  [34].  A  general  agreement  to  this  trend  is  stated  by  Hanazaki,  presented  as 
Figure  8  in  that  work  [21].  However,  the  Hanazaki  work  suggests  the  decrease  in  drag  coefficient 
begins  to  occur  after  Fd  =  0.25. 

In  turbulent  stratified  flows  around  a  sphere,  the  wake  region  can  be  broken  down  into  three 
regions:  the  near  wake,  the  non-equilibrium  regime,  and  the  far  wake  (e.g.  Spedding  1997  [49]  ). 
The  near-wake  is  located  where  the  flow  exhibits  similar  characteristics  to  the  unstratified  case. 
After  the  near-wake  regime,  there  is  the  so-called  non-equilibrium  (NEQ)  region,  where  buoyancy 
effects  become  apparent  as  the  flow  begins  to  exhibit  a  reduction  in  vertical  velocity  fluctuations. 
Finally,  there  is  the  far- wake  quasi-2D  regime  (Q2D),  wrhere  the  flow  characteristics  become  quasi 
two-dimensional  as  ’’pancake”  eddies  form.  The  near- wake  region  typically  exists  for  Nt  <  2;  the 
NEQ  region  for  2  <  Nt  <  50,  and  the  Q2D  regime  from  Nt  >  50  according  to  Spedding  1997  [49]. 

The  stratification  is  believed  to  have  little  influence  on  the  turbulence  of  the  near-wake.  The 
phenomenological  argument  is  made  that  in  the  very  near  wake  at  high  Reynolds  number,  the 
kinetic  energy  in  the  wake  is  high  compared  to  the  potential  energy  required  to  move  a  fluid  particle 
vertically  in  the  wake.  This  is  the  characteristic  of  the  near-wake  behavior.  The  work  by  Hopfinger 
et  al.  23]  gives  interesting  insight  into  the  behavior  of  the  flow  at  a  Reynolds  number  of  3000. 
The  Froude  number  parameter  denotes  a  possible  change  in  the  mechanism  for  generating  internal 
waves.  If  the  Froude  number  is  less  than  2,  then  the  wake  is  dominated  by  lee  waves.  If  the  Froude 
number  is  greater  than  2,  then  the  internal  wave  field  is  dominated  by  the  turbulent  wake  and 
internal  waves  produced  are  seemingly  random.  This  suggests  that  even  at  high  Reynolds  number 
flows,  the  stratification  has  the  ability  to  influence  the  flow  characteristics  even  when  stratification 
is  not  considered  severe. 

As  the  kinetic  energy  from  the  near-wake  decays,  the  available  surplus  in  converting  kinetic  to 
potential  energy  decreases  and  as  a  result,  the  growth  rate  of  the  wake’s  vertical  height  decreases 
so  that  the  wake  spreads  faster  in  the  horizontal  than  the  vertical.  This  is  the  essence  of  the  NEQ 
regime.  Schooley  &;  Stewart  1963  [41]  first  noticed  this  behavior,  is  also  noted  by  Pao  <k  Kao  1977 
[38]  and  reviewed  by  Lin  &  Pao  1979  [30].  The  NEQ  region  is  often  referred  to  as  the  ”  collapse” 
of  the  wake.  This  phraseology  is  not  globally  accepted  as  an  accurate  description  of  the  physical 
behavior  in  this  region.  Once  the  vertical  velocity  perturbations  are  negligible  when  compared  to 
the  horizontal  ones,  the  wake  has  evolved  into  the  far- wake  region  according  to  Bonnier  et  al.  1998 
[5]. 

Recent  interest  in  continuously  stratified  flows  around  submersed  bodies  is  partially  focused  on 
turbulent  far-wakes  and  whether  the  far-wrake  exhibits  any  characteristics  dependent  on  the  shape 
of  the  submersed  body.  In  2004,  Meunier  and  Spedding  speculated  that  the  momentum  induced  by 
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a  towed  body,  and  not  the  particular  shape,  is  important  for  the  characteristics  exhibited  by  the  far 
wake  [35].  In  2006,  Meunier  and  Spedding  generalized  that  result  to  include  the  effective  momentum 
diameter  hypothesis  to  encompass  self-propelled  submersed  bodies  [36].  Thus,  all  submersed  bodies 
traveling  through  a  stratified  fluid,  whether  towed  or  self-propelled,  may  be  appropriately  scaled 
by  an  effective  momentum  diameter  [36].  As  a  result,  a  sphere  traveling  horizontally  and  steadily 
through  a  continuously  stratified  fluid  is,  by  geometry,  its  own  effective  diameter  and  the  canon  of 
the  sphere  is  confirmed. 

Recent  computational  investigations  have  been  typically  focused  oil  effects  of  stratification  in 
the  far  wake  behind  the  sphere.  Currently,  far-wake  numerical  simulations  require  an  initialization 
procedure  to  create  an  initial  turbulent  flow-field  to  begin  integration  of  the  governing  equations 
in  time  ([16]  [15]).  Dommermuth  et  al.  2002  [16]  initialized  far  wake  simulations  based  in-part 
on  the  results  of  Bevilaqua  &:  Lykoudis  1978  [3].  The  data  available  is  a  Gaussian  distribution 
of  the  average  turbulent  kinetic  energy  [3].  Diamessis  et  al.  2005  also  followed  this  initialization 
procedure  [15],  and  recent  researchers  have  also  used  a  similar  scheme  [7].  It  should  be  noted  that 
results  in  [3]  are  not  collected  from  stratified  flow  experiments. 

None  of  the  aforementioned  simulations  have  included  the  sphere  explicitly  within  the  domain. 
This  is  because  the  computational  cost  to  both  model  the  proper  fluid  mechanics  on  the  sphere  and 
reproduce  the  far-wake  is  prohibitively  expensive  with  current  computational  resources.  Through 
non-dimensionalization  and  algebraic  relationships,  the  regimes  become  separated  by  N  *  t,  or 
equivalently  Fd~l  x/D.  In  a  study  where  Fd  —  5,  for  example,  the  near-wake  regime  may  not  end 
until  roughly  20  body  lengths  downstream  of  the  body,  and  the  Q2D  regime  does  not  begin  until 
250  body  lengths  downstream.  By  comparison,  the  boundary  layer  at  the  stagnation  point  of  a 
sphere  at  a  Reynolds  number  of  10,000  is  roughly  one-hundredth  of  a  body  length  thin. 

Brucker  &;  Sarkar  2010  [7]  suggest  that  there  still  is  ”no  dataset  which  characterizes  the  full 
state  of  the  turbulence  in  the  near  wake  region.”  There  is  then  an  opportunity  to  explicitly  account 
for  the  sphere  within  the  computational  domain  and  provide  detailed  information  about  turbulent 
fluctuations  as  well  as  the  density  field  in  the  near-wake.  The  density  perturbation  field  in  particular 
is  challenging  to  collect  experimentally.  Bonneton  et  al.  1996  [4]  explore  the  behavior  of  density 
stratification  in  the  wake  of  the  sphere  and  show  a  strong  correlation  between  the  density  and 
the  vertical  velocity  perturbations.  There  is  also  work  by  Bonnier  et  al.  2000  [6]  collecting  time 
dependent  density  field  information  in  the  far  wake,  but  that  work  is  limited  to  measurement  via 
conductivity  probe.  Many  early  experiments  rely  on  shadowgraph  or  dye  dispersion  techniques 
(e.g.  [31])  not  able  to  detect  perturbation  velocities,  and  other  optical  techniques  have  difficulties 
in  collecting  information  about  the  three  dimensional  density  field. 

According  to  the  literature  (e.g. [27]  [15]  [7]),  the  only  known  simulation  to  explicitly  account  for 
the  sphere  within  the  computational  domain  while  modeling  the  physical  problem  under  discussion 
is  that  of  Hanazaki  in  1988  [21].  Whilst  the  Froude  number  range  is  considerable  in  Hanazakfis 
1988  work,  the  chosen  Reynolds  number  and  solution  method  negates  any  unsteady  effects,  let 
alone  the  required  turbulent  fields  of  the  near- wake  that  may  be  of  interest  to  far-wake  modelers. 
The  proposed  investigation  fits  an  unfilled  niche  between  existing  experimental  work  done  in  the 
near  wake  and  Q2D  regimes  and  the  investigations  of  computational  modelers  ill  the  far  wake. 
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Chapter  2 

Numerical  Method 


2.1  Numerical  Method 


The  numerical  model  is  solved  on  a  discretized  grid  that  is  fitted  to  a  fixed,  structured  body-fitted 
coordinate  system.  The  Cartesian  coordinates  xx  are  transformed  into  generalized  coordinates  & 
but  the  flow  components  ut  and  density  perturbation  p  in  the  governing  equations  remain  in  the 
Cartesian  coordinates.  This  is  the  so-called  partial  transformation  formulation  of  the  governing 
equations  (e.g.  Const antinescu  1998  [12]).  The  effective  result  of  this  transformation  is  that 
the  physical  grid  is  transformed  into  a  computational  grid  by  the  transformation  relations  such 
that  the  Cartesian  coordinate  derivatives  are  represented  in  the  transformed  coordinates:  x  =< 
x(£,7/,C),  C)>  C)  This  formulation  has  been  used  successfully  in  several  studies  of 

varying  geometric  configuration  [25],  [12],  [11],  [13],  14]. 

The  governing  equations  must  again  be  recast  to  reflect  the  new  body-fitted  coordinate  system. 
A  review  of  these  transformations  can  be  found  in  Appendix  A  of  Sortiropolous  1991  [45]  or  in 
Chapter  2  of  Liseikin  1999  [33]. 

The  Jacobian  of  the  general  transformation  between  coordinate  systems  is  denoted: 

7  _ 

Jij  ~  ae 

The  determinant  of  the  Jacobian  is  simply  denoted  as  J.  Here,  xt  represents  the  Cartesian  co¬ 
ordinates  and  represents  the  coordinates  of  the  transformed  computational  grid.  Equations 
(1.8)-(1.1Q)  are  transformed: 
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(2.2) 

(2.3) 

(2.4) 


9 


and  Vi  is  the  contravariant  velocity.  Also,  in  Eq.  (2.2),  the  Cartesian  unit  vector  has  been  dropped 
from  the  explicit  statement  of  the  equation  and  is  presumed  to  be  understood. 


2.1.1  Computational  Domain 

The  numerical  model  is  discretized  using  finite  difference  schemes,  generally  of  2nd  or  5th  order 
for  the  momentum  equations  and  1st  order  for  density.  Finite  difference  calculations  are  performed 
on  the  computational  grid.  For  this  study,  the  physical  grid  is  constructed  in  spherical  coordinates 
in  an  0-0  topology,  where  there  are  known  relations  between  <  r,  6,<p  >  and  x  such  that  there  are 
also  known  relations  between  x  and  £,  similar  to  what  is  done  in  [25] : 


r  = 


Rma 


D\  /tanh(fci  *u^T-1) 

2/1  tanh  fci 


(2.5) 


where  subsequently: 


0  =  7T- 


0  =  2tt 


g-1 

(max  —  1 

T)-  1 


V max  2 


x  =  r  cos  0,  y  —  r  sin  0  cos  0,  z  —  r  sin  0  sin  <p 


(2.6) 

(2.7) 

(2.8) 


Rmax  is  chosen  to  approximate  a  far-field  condition  in  the  simulations.  Figure  4.2  gives  a  graphical 
representation  of  a  typical  physical  grid  used  in  this  study.  A  sketch  of  the  body  fitted  coordinate 
system  as  related  to  the  sphere  is  found  in  Figure  4.3. 

The  computational  grid  is  then  constructed  from  the  known  relations  of  Eqs.  (2.5)-(2.8).  A 
visualization  of  the  computational  grid  is  in  Figure  4.4.  There  is  a  periodic  overlap  of  one  grid  point 
in  the  azimuthal  direction  such  that  rji  =  v/max- 1  and  rfr  =  rjmax-  Additionally,  a  grid  singularity  is 
created  by  the  choice  of  coordinates  in  the  physical  domain,  located  at  0  =  {0  =  £ min ,  n  =  £max}. 


2.1.2  Numerical  Solver 


The  governing  equations  are  integrated  by  an  added  transient  term  in  pseudo-time  such  that 
when  the  solution  exhibits  a  steady  state  in  pseudo-time,  the  governing  equations  have  also  con¬ 
verged  in  real  time.  This  is  the  main  concept  of  a  so-called  dual-time  stepping  scheme,  and  the 
implementation  by  Arnone  et  al\  1]  is  followed. 

Shorthand  notation  for  specific  operators  within  the  governing  equation  is  introduced: 


CONV(  )  =  Vj 


™£SS<  >  -  §tw 


v,sci  >  =  7S 


and  the  explicit  discretizations  of  the  continuous  operators  CONV(  ),  PRESS(  ),  and  V I  SC  (  ) 
are  found  in  Constantinescu  1998  [12]. 


10 


Flow  quantities  are  linearized  for  implicit  treatment  of  flow  variables  in  pseudo-time  integration. 


Qm+l  =  Qm  +  AQ 


(2.9) 


Then  the  momentum  equation  (2.2)  becomes: 


_ 1/ )i  ^ _ i/*^ 

‘  ■+  -it-1  +  CONVm(uim+l) 


At 


At 


/m+l 

+  PRESS(p'm+l)  -  VISC(uim+1)  +  ^rp-h 3  =  0 


(2.10) 


”m”  indicates  the  m-th  discrete  point  in  pseudo-time  integration  and  ”n”  indicates  the  previous 
level  in  the  discretized  physical  time  integration.  Each  iteration  begins  with  m  =  n,  and  as 
limAQ— o  Qm+l  =  Qm,  implying  Qm+1  =  Qn+1. 

The  momentum  equation  is  split  in  psuedo-time  to  solve  for  both  the  momentum  and  pressure 
variables.  The  details  are  available  in  Constantinescu  1998  [12].  Once  the  momentum  and  pressure 
field  are  projected  to  level  ”m-Fl’\  the  additional  equations  for  the  turbulence  model  and  the 
density  perturbation  are  advanced  in  pseudo-time.  The  superscript  in  Eq.  (2.10)  for  the  density 
variable  is  denoted  level  ”m”  for  this  reason. 

The  solver  uses  a  non-staggered  formulation  that  stores  both  pressure  and  velocity  at  the  nodal 
grid  points.  Without  a  proper  formulation,  this  can  create  the  now-famous  ’’checkerboarding” 
problem.  To  counter  this  potential  issue,  the  solver  uses  the  method  of  Sortiropolous  and  Patel 
[46],  which  adds  an  artificial  dissipation  term  to  the  pressure  equation.  Values  of  the  pressure 
dissipation  coefficient  used  in  this  study  were  usually  0.05-0.1.  The  details  of  its  implementation 
within  the  numerical  integration  scheme  can  be  found  in  Constantinescu  1998  [12]. 


2.1.3  Turbulence  Modeling 

As  discussed  in  Section  1.2.2  the  Reynolds  number  ranges  investigated  in  this  proposal  are  antic¬ 
ipated  to  exhibit  turbulent  wake  behavior.  Because  of  the  finite  discretized  nature  of  the  numerical 
experiment  and  present  limits  on  computational  resources,  the  turbulent  motions  must  be  approx¬ 
imately  modeled.  Turbulence  modeling  is  a  still- active  area  of  research,  and  discussion  will  be  kept 
to  the  type  of  modeling  employed  in  this  numerical  experiment. 

Reynolds  Averaged  Navier  Stokes  (RANS)  formulations  are  an  extensively  used  method  to 
model  turbulent  flows.  They  arc  particularly  suited  for  flows  where  only  the  mean  flow  properties 
are  desired  (c.g.  Wilcox  pg.  30  [55]).  RANS  formulations  are  based  on  the  idea  of  averaging  flow 
field  properties,  and  they  tend  to  ”  average-out”  the  temporal  features  ill  the  flow  as  the  number 
of  samples  or  length  of  averaging  time  increases.  RANS  formulations  are  particularly  useful  in 
exploring  the  behavior  of  separated  boundary  layers;  assuming  the  separation  point  exhibits  a 
quasi-steady  nature.  RANS  formulations  are  unsuitable  if  the  time-dependent,  spatially  developing 
structures  of  the  flow  away  from  the  boundary  are  of  interest,  especially  in  the  axisymmetric  case 
of  the  sphere  [13]. 

An  alternative  approach  that  attempts  to  represent  the  aforementioned  large-scale  motions  is 
referred  to  as  Large  Eddy  Simulation  (LES).  Large  eddy  simulation  attempts  to  spatially  and 
temporally  resolve  large-scale  motions  in  space  and  time  but  to  model  any  motion  that  is  spatially 
unresolvable.  These  spatially  unresolvable  motions  are  commonly  referred  to  as  the  sub-grid  scale 
motions.  A  major  drawback  of  LES  occurs  in  modeling  flows  in  the  presence  of  a  wall,  and 
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consequently  boundary  layers.  This  is  because  the  scales  that  exist  within  boundary  layers  are 
very  small  compared  to  the  other  large-scale  motions  in  the  flow.  As  such,  the  spatial  resolution 
required  in  pure  LES  is  very  fine,  and  computational  costs  can  be  comparable  to  Direct  Numerical 
Simulation;  which  is  very  costly.  Thus  LES  flows  near  boundaries  can  be  cost-prohibitive  and 
current  research  in  LES  simulations  deals  with  wrall- modeling  close  to  boundaries. 

Instead  of  RANS  or  LES  alone,  Detached  Eddy  Simulation  (DES)  is  used  in  this  investigation. 
The  essence  of  DES  is  that  it  uses  RANS  as  a  wall-model  within  boundary  layers  and  an  LES 
sub-grid  scale  type  model  for  resolution  of  the  flow  away  from  walls.  DES  has,  on  its  own,  a  variety 
of  formulations  (e.g.  Spalart  et  al.  1997  [48]  or  Strelets  2001  [50])  and  is  not  specific  to  a  particular 
turbulence  modeling  formulation;  only  that  it  behaves  in  the  prescribed  RANS/LES  manner. 

The  current  investigation  uses  the  Spalart-Allmaras  DES  (SA-DES)  formulation  of  Constanti- 
nescu  2004  [14].  Near  solid  boundaries,  where  separation  occurs,  SA-DES  uses  the  one-equation 
Spalart-Allmaras  RANS  formulation  [47]  to  predict  boundary  layer  behavior;  including  separation. 
Outside  of  the  boundary  layer  region,  the  model  automatically  switches  to  an  eddy  resolving  LES- 
type  mode.  The  LES-type  model  behaves  like  a  Smagorinsky  model,  which  in  itself  attempts  to 
represent  the  Kolmogorov  spectrum  for  sub-grid  scales  as  determined  by  Lilly  in  1967  [29].  With 
this  formulation,  computational  costs  can  remain  closer  to  LES  even  with  the  sphere  present  in 
the  simulation.  It  should  also  be  noted  that  the  SA-DES  eddy  viscosity  may  be  used  for  the  den¬ 
sity  perturbation’s  eddy  viscosity  since  the  both  the  molecular  Schmidt  number  and  the  turbulent 
Schmidt  number  are  taken  to  be  unity  in  the  turbulent  cases  as  is  alluded  to  in  Section  1.1.2. 

The  addition  of  the  DES  turbulence  model  to  the  equations  of  motion,  Eq.  (2.10),  modifies  the 
viscous  operator  to  be  time  dependent;  i.e.  VISC(  )  — ►  VISCm(  )  in  the  following  way: 


where  =  v\n /{UD),  where  vt  is  the  eddy  viscosity. 

Consequently,  Eq.  (2.10)  stays  effectively  unmodified  to  the  user,  and  the  turbulence  model  is 
turned  off  as  the  term  — >  0.  The  SA-DES  formulation  for  the  eddy  viscosity  is  reproduced 

from  Nikitin  et  al.  2000  [37]  as  follows: 


where: 


%  +  V1§&  =  rbl^  +  “  [(V  •  (^  +  I5)  V  ^)  +  0.2  (V0)2] 


Cwlfw  \  ^ 


Pi2 


Vt  =  uf„l,  /^l  = 


v 

X  =  - 


X3  +  4  1  ’  v 

S  =  S  +  (t> / kcP)  ,  /„2  =  1  -  ,\/(l  +  xM 


fw  =  9 


_L+ik16 

g6  +  4.3 


,  g  =  r  +  c.u2  (r6  -  r)  ,  f  = 


SK2d 


d  -  min  ( dw ,  Cdes&) 


(2.11) 
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and  S  is  the  magnitude  of  the  vorticity;  dw  is  the  physical  distance  to  the  wall;  A  is  the  local  grid 
spacing;  cb\  =  .1355,  a  =  2/3,  cb 2  =  0.622,  k  =  .41,  cwX  =  chx/K2  +  (1  +  c/?2)/cr,  =  0.3,  c^3  =  2, 

and  c„i  =7.1  are  modeling  constants.  The  value  of  Cdes  is  generally  set  at  0.65  following  the 
work  done  by  Shur  et  al.  [44]  for  homogeneous  turbulence. 


2.1.4  Boundary  Conditions 

Momentum,  Turbulence  Model,  and  Density  B.C.’s 

On  the  surface  of  the  sphere,  Rmin  =  Cmm,  Dirichlet  boundary  conditions  are  enforced  for  momen¬ 
tum  and  the  turbulence  model  variable,  whilst  density  perturbations  use  the  zero-normal  gradient 
boundary  condition  to  enforce  incompressibility  within  the  diffusive  density  equation: 


-  dp 

o,  £  =  0,  -^  =  0 

on 


(2.12) 


At  the  domain  boundary  far  from  the  sphere,  mixed  boundary  conditions  are  used  due  to  the 
geometry  of  the  physical  grid.  The  first  condition  is  referred  to  as  an  ” inflow”  boundary  and  has 
a  Dirichlet  condition  for  the  velocity,  the  turbulent  variable,  and  any  density  perturbations  for  the 
surface  along  {£min  =  0  <  9  <  .557T,  Cm ou:  =  Rmax}‘- 


u  =  cl,  v  =  0,  p  —  0 


(2.13) 


If  the  inflow  boundary  is  far  upstream  from  the  sphere,  it  is  considered  sufficient  to  use  the  Dirichlet 
conditions  for  low  Reynolds  number,  stratified  cases,  according  to  Hanzaki  1988[21].  Hanazaki  used 
an  outer  domain  distance  of  20  sphere  diameters  [21].  In  the  current  study,  the  domain  boundary  is 
set  14.5  sphere  diameters  away  from  the  surface  of  the  sphere.  This  distance  was  found  sufficient  for 
low  Reynolds  simulations  in  [25]  and  was  also  found  sufficient  for  turbulent,  high  Reynolds  number 
(Re  >  105)  simulations  [14].  Both  [25]  and  [14]  were  simulations  of  a  sphere  traveling  steadily 
through  a  fluid  of  a  uniform  density. 

The  second  condition  on  the  domain  boundary  awray  from  the  sphere  is  considered  an  ” outflow” 
boundary  condition.  Outflow  boundary  conditions  rely  on  information  from  inside  the  compu¬ 
tational  domain.  The  outflow  boundary  is  contained  on  the  surface  where  {.557T  <  6  <  tt  = 
CmajM  Cm  ax  =  Rm  ax  } : 


D  14, 

Dt 


+  j£pfa-0’ 


Dp_ 

Dt 


u3  =  0 


(2.14) 


For  momentum,  these  conditions  allow  a  fluid  particle  to  accelerate  in  the  vertical  duo  to  buoyancy 
forcing  but  without  accelerating  in  the  horizontal.  Incompressibility  is  applied  via  the  density 
equation;  leaving  u3  in  the  formulation  of  the  boundary  condition  to  represent  the  transport  of  the 
background  fluid  density  (e.g.  Kundu  pg.  249  [26]). 

On  the  singularity  lines  present  at  9  =  {0,7r},  all  variables  for  momentum,  the  turbulence 
model,  and  density  are  extrapolated  from  the  interior  and  averaged  over  the  repeated  nodes.  Along 
the  periodic  overlap  in  the  ?/  direction,  periodicity  is  applied  to  each  respective  flow  variable. 


Pressure  B.C.’s 

The  pressure  boundary  condition  is  controlled  by  using  the  normal  component  of  the  momen¬ 
tum  equation  Eq.  (2.2)  on  the  computational  boundaries  located  at  Rmin  and  Rinax  to  enforce 
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conservation  of  mass  [19]: 


*C  (2.15) 

This  provides  an  explicit  equation  for  the  pressure  variable  when  put  into  discretized  form  along 
both  the  sphere  surface  and  the  far  field  domain  boundary. 

Along  the  singularity  lines  present  at  0  =  {0,7r},  the  pressure  is  extrapolated  from  the  interior; 
subject  to  the  conservation  of  mass  constraints  and  similar  to  what  is  recommended  in  Johnson 
1999  [25].  Along  the  periodic  overlap  in  the  q  direction,  periodicity  is  applied  to  the  pressure  field 
as  well. 


V(P)  =  ~^~  CONV(u)  +  VISC(u)  -  -jL 


2.1.5  Computational  Outflow  Boundary  Treatments 

The  physical  problem  being  modeled  exists  within  an  infinite  domain  and  may  contain  internal 
waves  that  propagate  from  within  the  turbulent  wake  structure  (see  Section  1.2.2).  An  infinite 
domain  would  allow  any  internal  waves  to  propagate  out  towards  infinity,  which  is  the  expected 
behavior  (e.g.  Lighthill,  Section  4.9  [28]).  This  is  known  as  the  radiation  condition. 

Introducing  an  artificial  boundary  at  the  edge  of  the  physical  domain  complicates  computational 
simulations  containing  waves  because  the  radiation  condition  may  not  be  appropriately  met.  Tur¬ 
bulent  wakes  produce  waves  close  to  the  computational  boundary.  The  exact  correction  required  to 
ensure  the  radiation  condition  to  be  completely  satisfied  is  as-yet  elusive,  and  several  approaches 
are  actively  under  investigation  (e.g.  Givoli  2004  [18]). 

Viscous  Damping 

It  is  desirable  to  find  a  method  capable  of  reducing  these  potentially  reflected  waves  while  at 
the  same  time  maintaining  a  certain  simplicity  in  its  implementation.  A  tempting  methodology 
is  to  introduce  some  artificial  viscous  damping  to  the  equations  of  motion,  Eq.  (2.2),  near  the 
region  of  the  boundary,  but  far  enough  from  the  region  of  interest  so  as  not  to  adversely  affect  the 
investigation’s  observations: 

NS(q)  =  -v(x)q  (2.16) 

NS(q)  =  fi(x )0  (2.17) 

where  NS(q)  represents  the  original  Navier-Stokes  bracketed  term  in  Eq.  (2.2)  operating  on  any 
desired  quantity  q  that  is  to  be  damped.  i/(x)  and  /i(x)  are  the  damping  coefficients  which  are 
allowed  to  vary  spatially.  Israeli  and  Orszag  1981  [24]  show  that  the  damping  term  must  vary 
smoothly  in  space  to  avoid  reflections  off  of  the  damper  interface  itself,  and  they  show  that  Eq. 
(2.16)  is  superior  in  its  wave  damping  properties.  In  flows  with  turbulent  quantities  where  fluctua¬ 
tions  may  be  centered  about  some  mean  profile  value,  Eq.  (2.16)  may  be  extended  to  what  is  used 
in  Brucker  2010  [7]: 

NS{q)  =  -v{x){Q-q)  (2.18) 

where  Q  represents  the  characteristic  mean  quantity  of  the  flow  variable  q  that  is  to  be  attained 
within  the  damping  region. 
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Damping  with  Zonal  DES 

An  original  unpublished  method  that  may  be  viable  is  to  use  the  eddy  viscosity  of  the  turbu¬ 
lence  model  itself,  rather  than  an  arbitrary  addition  to  the  equations  of  motion,  to  damp  out  any 
oscillations  towards  or  near  the  boundary.  The  idea  is  to  use  the  ’’switch”  inherent  to  the  DES 
formulation  Eq.  (2.11)  to  explicitly  turn  off  the  LES  mode  towards  the  outflow  boundary  and 
return  the  turbulence  modeling  to  RAXS;  explicitly: 

d  =  dw  for  x  e  M(xdamp)  (2.19) 

where  dw  is  the  physical  distance  to  the  wall  as  in  the  original  SA-RAXS  model  [47]  and  Q(xdamp) 
is  the  domain  in  which  damping  is  desired.  This  w’ould  create  a  RAXS  simulation  at  the  outflow 
and  would  maintain  the  DES  behavior  desired  in  the  separation  and  wake  region. 

An  additional  possibility  of  using  the  DES  formulation  as  a  damper  is  to  increase  the  value  of 
Coes  found  in  Section  2.1.3  for  areas  near  the  computational  outflow  boundary.  Constaiitinescu 
and  Squires  2004  [11]  show  that  by  increasing  the  value  of  Cdes  the  smaller  scales  in  the  energy 
spectra  decay  as  Coes  is  increased.  Conceptually,  for  high  enough  vales  of  Coes ,  damping  acts 
similar  to  RAXS. 
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Chapter  3 


Preliminary  Results 

3.1  Laminar  Regime 

3.1.1  Re=200  DNS  Results 

In  order  to  ensure  the  code  produces  viable  results  for  stratified  flow  cases,  it  is  important  to 
compare  code  results  with  existing  data  sets.  Since  Hanazaki  1988  [21]  is  related  to  the  proposed 
work,  it  is  a  natural  benchmark  to  compare  the  current  numerical  code  against. 

Perhaps  a  natural  starting  point  is  in  comparison  of  the  drag  coefficients.  Results  in  Figure  4.8 
are  in  excellent  agreement  with  those  of  Hanazaki,  except  for  the  low  Fd  —  0.125  case.  This  is  an 
interesting  discrepancy  between  the  two  works  as  they  tend  toward  agreement  in  all  other  cases  at 
Re  =  200. 

Figure  4.7,  which  is  a  measure  of  the  change  in  stream- wise  drag  coefficient  as  function  of  the 
Froude  number  at  a  Reynolds  number  of  200.  Here,  after  a  value  of  Fd  =  0.25  (1  /Fd  =  4),  the 
data  of  Hanazaki  sees  a  sharp  decline  in  the  drag  coefficient  as  the  affect  of  stratification  on  the 
surrounding  flow  becomes  more  severe,  whereas  the  value  of  the  current  work  trends  to  capture 
experimental  trends  more  closely.  The  experimental  data  of  both  Mason  and  Lofquist  and  Purtell 
show  an  increasing  trend  until  Fd  ~  0.15  (1  /Fd  =  7);  followed  by  a  gradual  decline  in  ACd  as 
1/F  approaches  a  value  of  10. 

Hanazaki  suggests  that  the  discrepancy  between  the  computational  results  and  prior  experi¬ 
mental  data  is  due  to  Reynolds  number  difference.  The  experimental  data  of  Lofquist  &  Purtell  84 
[34]  show  that  experiments  were  conducted  at  a  Reynolds  number  of  approximately  220  for  the  low 
Froude  number  cases.  This  is  not  a  far  departure  from  the  numerical  simulations  at  the  Reynolds 
number  of  200  reported  by  Hanazaki.  Hanazki  does  not  report  on  vortex  shedding  present  at  these 
low  Reynolds  numbers  [21],  an  observation  also  noted  by  Chomaz  et  al.  1993  [9]. 

The  simulations  presented  in  this  investigation  do  capture  the  vortex  shedding  behavior  regime. 
Analysis  of  the  drag  coefficient  curves  presented  in  Figure  4.9  shows  that  there  is  a  strong  lateral  (y- 
plane)  periodic  force  indicative  of  vortex  shedding.  Spectral  analysis  of  the  lateral  drag  coefficient 
signal  shows  that  there  is  a  Strouhal  number  of  0.18  for  the  Fd=0.125  case.  Strouhal  numbers  for 
this  case’s  drag  signals  can  be  seen  in  Table  3.1. 

The  shedding  of  the  vertical  component  of  vorticity,  u?Zl  over  the  course  of  1  shedding  cycle  is 
presented  in  Figure  4.11.  Qualitatively,  it  is  apparent  that  the  sign  of  urz  on  either  side  of  the  x-axis 
alternates  every  quarter  period  in  the  wake  region. 

In  homogeneous  fluid  flows  around  a  sphere,  drag  coefficients  may  also  depend  on  the  location 
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St 

Power 

x: 

y- 

z: 

{0.035,0.355} 

0.18 

{28, 16.23} 
426.5 

No  significant  peaks. 

Table  3.1:  Drag  coefficient  Strouhal  numbers:  Re=200,  Fd=0.125 


of  the  separation  point/line  on  the  sphere.  To  locate  the  separation  point  on  an  azimuthal  plane 
of  the  sphere,  the  skin  friction  coefficient  is  used  as  an  indicator  for  separation.  The  skin  friction 
coefficients  for  a  range  of  Fd  are  shown  in  Figure  4.12. 

As  mentioned  in  Section  1.2.1,  a  zero  skin  friction  coefficient  may  not  be  the  best  absolute 
indicator  for  separated  flow.  A  qualitative  idea  of  what  is  physically  occurring  in  these  steady  flows 
can  be  derived  by  looking  at  the  streamlines  close  to  the  sphere  as  in  Figure  4.14.  The  necessity  of 
this  additional  information  may  be  seen  in  in  Figure  4.12  (c),  where  there  is  an  apparent  crossing  of 
the  zero  value  for  polar  values  of  roughly  145°  and  160°  for  the  vertical  azimuthal  plane.  Inspection 
of  Figure  4.14  (c)  shows  that  particles  in  the  flow  eject  away  from  the  body  of  the  sphere,  and  indeed 
the  separation  point  is  located  at  145°  where  this  first  occurs. 

In  Figure  4.13,  comparison  is  made  of  the  vertical  and  horizontal  separation  points’  location 
as  compared  to  the  data  of  Hanazaki  1988.  It  becomes  apparent  from  this  figure  that  there  are 
differences  between  the  two  cases. 

As  suggested  in  Section  1.2.2,  one  of  the  main  advantages  of  computational  investigations  into 
stratified  flows  is  the  relative  ease  which  it  is  possible  to  obtain  instantaneous  density  profiles  of 
the  entire  flow  field.  In  Hanazaki  1988,  isopycnals  are  presented  and  that  trend  is  followed  here. 

An  advantage  of  the  linear  density  stratification  is  that  with  non-dimonsionalization,  the  actual 
density  gradient  term  ’’drops  out”  of  the  governing  equations.  Thus,  any  linear  gradient  may  be 
imposed  on  the  density  perturbation  field  to  represent  the  more  physically  intuitive  isopycnals; 
the  numerical  values  of  which  are  then  arbitrary.  Isopycnals  from  the  simulations  are  presented  in 
Figure  4.15.  From  this  figure,  the  flow  slowly  becomes  more  two  dimensional  as  the  Froude  number 
decreases;  signifying  the  tendency  of  the  flow  to  begin  its  pseudo  two  dimensional  layered  behavior. 

Finally,  contours  of  vertical  components  of  velocity,  w,  are  presented  in  Figure  4.16.  These 
contours  for  this  low  Reynolds  number  case  serve  as  contrast  contours  to  identify  lines  of  constant 
phase.  Lines  of  constant  phase  are  located  where  the  vertical  component  of  velocity  is  zero  according 
to  linear  theory  as  u)  —  N  *  sin(0)  [28].  Here,  0  is  the  angle  of  the  motion  of  the  fluid  particle  with 
respect  to  the  horizontal;  making  lines  of  w  =  0  equivalent  to  lines  of  constant  phase. 

Perhaps  the  most  striking  feature  in  comparing  Figure  4.16  (a)-(d)  is  that  the  lines  of  constant 
phase  appear  to  have  similar  patterns  in  the  free  stream  portion  of  the  flow,  but  different  spacing 
depending  on  Froude  number.  Figure  4.16  hints  that  reported  values  for  flow  behind  the  sphere  in 
stratified  flow  simulations  should  perhaps  be  given  in  terms  of  possibly  ^  versus  fy. 

3.1.2  Re=1000  Preliminary  Results 

Comparisons  are  initially  made  with  the  laminar,  unsteady  unstratified  case  at  a  Re=300  of 
Johnson  [25]  to  verify  the  unsteady  behavior  of  the  numerical  solver.  The  resulting  stream-wise 
drag  coefficient  in  the  prepared  case  is  Cdx  =  .654±.004,  Cdy  —  —  .054±.015,  and  Cdz  =  —,035±.01. 
As  is  shown  in  Table  3.2,  the  average  drag  coefficient  in  the  stream- wise  direction  and  the  drag 
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coefficient  on  a  rotated  plane  of  symmetry  (the  lateral  plane,  denoted  by  the  L  subscript)  compare 
favorably  to  Johnson’s  results. 


Orr  ’10 

Johnson  '99 

Cdz 

0.654  ±  0.004 

0.656  ±  .0035 

Cdl 

0.064  ±  0.019 

0.069  ±  .0016 

Stx 

0.13 

0.137 

StL 

0.13 

0.137 

Table  3.2:  Drag  Coefficient  and  Strouhal  Number  at  Re=300,  Fd=oc 


Orr  10 

Lee  ’00 

cDl 

0.555  ±  0.04 

0.545 

CDy 

-0.022  ±  0.08 

N/R. 

Cd, 

-0.017  ±  0.155 

N/R 

Stx 

{0.04  to  0.05,  0.135  to  0.15} 

{0.043,  0.164} 

Table  3.3:  Drag  Coefficient  and  Strouhal  Number  at  Re=500,  Fd=so.  N/R:  not  reported 

Further  comparison  to  the  work  of  Lee  2000  [27]  is  made  to  compare  the  unsteady  behavior  at 
a  Reynolds  number  of  500.  The  time  series  drag  coefficient  data  is  again  compared  in  Table  3.3. 
Results  appear  to  compare  favorably  against  this  data  set.  Strouhal  numbers  associated  with  both 
the  unsteady  motion  of  the  separation  point  and  the  unsteady  vortex  shedding  produce  similar 
temporal  behavior;  within  a  10%  error  of  what  is  produced  by  Lee.  There  are  inherent  differences 
between  the  numerical  method  of  Lee  2000  and  the  current  study,  and  it  is  possible  that  either 
grid  resolution,  or  the  discretization  scheme,  or  even  the  solution  method  may  account  for  this 
discrepancy.  Of  interest  is  that  the  ratio  of  peak  Strouhal  numbers  between  Lee  and  the  current 
study  is  within  a  few  percent;  suggesting  that  whatever  the  discrepancy  in  exact  value  of  Strouhal 
numbers,  the  relations  between  the  various  physical  mechanisms  of  the  simulations  seem  to  be  in 
agreement. 

After  the  above  comparisons  to  check  unsteady  behavior,  the  Reynolds  number  is  further  in¬ 
creased  to  1000.  The  general  approach  for  Re=200  is  followed  as  a  guide  for  the  analysis  of  Re=1000. 
At  this  Reynolds  number,  there  are  four  distinct  regimes  of  the  wake  reported  by  Lin  et  ai  [32]: 
the  lee  wave  instability,  the  non-axisymmetric  vortex,  the  symmetric  vortex  shedding,  and  the  non- 
symmetric  vortex  shedding.  Data  is  presented  at  Froude  numbers  of  2,  0.5,  and  0.25.  These  points 
correspond  to  non-symmetric  vortex  shedding,  non-axisymmetric  vortex,  and  lee  wave  instability 
regions,  respectively.  The  data  for  the  expected  symmetric  vortex  shedding  region  is  not  available 
at  this  time. 

The  relationship  for  the  drag  coefficients  as  seen  in  Figure  4.17  remain  relatively  unchanged  from 
the  Re=200  cases  presented  in  the  Figure  4.7  counterparts.  It  will  be  of  interested  to  examine  the 
Fd=0.125  case  to  see  if  the  trend  at  Re=1000  is  continued  as  is  seen  at  Re=200  for  the  stream- wise 
drag  coefficient. 

Inspection  of  the  Re=1000,  Fd=0.25  case  warrants  further  explanation  as  unsteady  effects  were 
observed  whereas  the  guide  of  Lin  et  ai  suggests  a  steady  lee  wave  formation  behind  the  sphere. 


18 


Evidence  of  unsteady  vortex  shedding  is  seen  by  inspection  of  the  drag  coefficient  time- series  in 
Figure  4.22  and  supported  by  the  strong  spectral  peak  contained  in  Figure  4.23.  Qualitative 
appearance  of  vortex  shedding  similar  to  the  behavior  in  Figure  4.11  has  been  confirmed  and  is  not 
reproduced  here. 

Inspection  of  the  pressure  coefficients  in  the  vertical  (z  >  0)  and  horizontal  ( y  <  0)  in  Figure 
4.18  suggest  that  while  general  relationships  between  the  two  planes  remain  roughly  the  same,  there 
are  a  few  quantitative  differences  between  the  Re=200  and  Re=1000  case  as  the  Froude  number 
decreases.  This  is  especially  clear  at  a  Froude  number  of  0.25,  where  an  internal  hydraulic  jump  (as 
defined  by  Hanazaki  1988  [21]  is  seen  for  the  Re=1000  ease,  whereas  it  is  not  seen  for  the  Re=200, 
Fd=0.25  case. 

By  comparing  Figure  4.12  and  Figure  4.19,  interesting  agreement  is  seen  in  the  behavior  of  the 
separation  point  at  a  Reynolds  number  of  1000  is  qualitatively  similar  to  the  Reynolds  number  of 
200  case  for  varying  Froude  numbers.  The  separation  point  is  notably  delayed  in  the  vertical  plane 
by  the  stratification.  The  separation  point  also  does  not  seem  to  move  as  rapidly  downwind  in  the 
horizontal  as  compared  to  the  vertical  planes,  just  as  in  the  Reynolds  number  200  case.  This  may 
be  inferred  from  Figure  4.20,  but  a  more  diverse  range  of  Froude  numbers  for  data  collection  would 
be  prudent.  Comparison  between  Figure  4.13  and  Figure  4.20  is  a  challenge  due  to  lack  of  available 
data  at  Re=1000.  It  does  appear  that  the  variance  of  vertical  and  horizontal  separation  point 
decreases  between  the  two  Reynolds  numbers,  which  may  owe  to  the  higher  energy  contained  in  the 
higher  Reynolds  number  flows.  The  data  point  at  Re=200,  Fd=0.5  for  the  horiztonal  separation 
has  been  verified  as  a  simulation  result,  but  the  reason  for  this  ’’odd”  data  point  location  is  not 
yet  known.  Separation  points  may  be  qualitatively  inferred  from  inspection  of  the  instantaneous 
velocity  vector  field  as  the  boundary  layer  point  at  separation  is  assumed  to  be  quasi-steady.  This 
is  available  in  Figure  4.21. 

The  largest  difference  between  the  Re=1000  cases  and  the  Re=200  cases  is  the  ability  to  observe 
unsteady  effects.  Time  series  signals  are  then  used  as  the  analysis  tool,  particularly  for  the  density 
perturbation  field.  A  probe  is  placed  at  an  Fd~l  x/D  ~  4  (Fd=2)  to  observe  the  unsteady  behavior 
of  the  flow  in  the  near-wake.  By  inspection  of  Figure  4.24  and  Figure  4.25  it  is  possible  to  see  that 
there  is  a  clear  correlation  between  the  vertical  velocity,  wf,  and  the  density  perturbation  field,  p \ 
Note  that  inspection  between  Figure  4.24  (a)  and  (b)  at  t*  ~  130  show  a  slight  difference  in  signal 
values,  indicating  they  are  not  perfectly  correlated,  but  the  two  signals  do  remain  highly  correlated. 
Figure  4.25  serves  as  an  additional  qualitative  indicator  of  this  and  clearly  w f  is  a  leading  order 
term  in  Equation  1.10,  but  analyses  of  the  phase  lag  of  the  two  signals  may  be  required  as  complete 
evidence  of  this  claim. 

Inspection  of  the  lines  of  constant  phase  in  Figure  4.26  indicates  that  the  variance  in  Reynolds 
number  does  not  appear  to  affect  the  internal  wave  properties  in  the  free  stream  as  is  also  seen  in 
Section  3.1.1. 
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3.2 


Turbulent  Regime 


3.2.1  Re=5000  SA-DES  Results 

Based  on  available  data,  this  investigation  is  also  performed  at  a  Reynolds  number  of  5000, 
Froude  number  of  2,  and  Schmidt  number  of  1  while  using  the  SA-DES  formulation  discussed  in 
Section  2.1.3.  All  discussion  within  this  section  uses  these  parameters.  Two  sets  of  data  are  used 
to  benchmark  the  behavior  of  the  simulation  at  this  Reynolds  number. 

First,  the  perturbation  velocities  of  the  simulations  axe  measured  against  the  current  computa¬ 
tional  standard  procedure  of  using  the  azimuthally  averaged  velocities  u\vf,w'.  Both  the  average 
and  perturbation  velocities  of  the  simulation  is  compared  against  the  Gaussian  profiles  of  Diamessis 
et  al.  [15].  All  data  is  averaged  in  the  azimuthal  direction,  </>,  then  spatially  averaged  along  the 
x-direction.  The  radial  perturbation  velocities  here  are  defined  as: 

u'(r)  =  ((u(x,t)  -  14(0)*)  z 
v’(r)  =  ([v(x,t))^x 
w'{r)  =  ((w{x,t))^ 

Ux(r)  = 

where  u,  v ,  w  represent  the  Cartesian  vector  components  of  the  velocity  vector  field  and  Ux(r) 
is  the  stream- wise  averaged  velocity  dependent  only  on  radial  distance  from  the  x-axis.  uf  is  also 
referred  to  as  the  stream-wise  perturbation  velocity. 

The  data  is  collected  at  the  same  location  that  Diamessis  et  al.  initialize  their  simulations  [15] 
at  Nt  =  3,  which  is  the  equivalent  downstream  distance  of  (1  /Fd)(x/D)  =  3,  taken  from  the 
center  of  the  sphere.  The  mean  streamwise  velocity  profile  is  compared  in  Figure  4.27.  The  mean 
average  velocity  compares  relatively  well  to  the  Gaussian  profile  fit.  This  appears  to  indicate  that 
the  SA-DES  model  does  a  relatively  good  job  of  predicting  the  mean  flow  properties  in  the  wake. 

Further  comparison  is  made  against  the  Gaussian  fitted  profile  of  [15],  where  the  perturbation 
velocities  u',  v\  wf  are  averaged  together.  This  is  the  assumption  made  for  the  velocity  perturbation 
initial  condition  field  by  Dommermuth  et  al.  [16]  and  Diamessis  et  al.  [15]  as  well  as  Bruker  & 
Sarkar  [7],  The  comparison  can  be  seen  in  Figure  4.29.  There  is  the  expected  drop  in  turbulent 
fluctuations  at  the  center  of  the  core,  indicating  a  relatively  lower  region  of  turbulent  production 
there.  This  trend  is  also  reflected  by  the  Gaussian  fit.  Again,  the  computational  result  appears  to 
compare  favorably  to  the  Gaussian  averaged  profile  currently  used  by  computational  modelers. 

The  individual  perturbation  velocities  are  also  investigated.  The  result  can  be  seen  in  Figure 
4.29.  As  seen,  the  perturbation  velocities  are  not  equally  contributing  to  the  average  profile  seen  in 
Figure  4.29.  The  components  of  the  perturbation  velocities  make  up  the  turbulent  kinetic  energy, 
i.e.  k  =  1/2  (u2  4*  v2  4-  w 2)^2,  and  thus,  the  question  is  raised  of  whether  the  assumption  of  the 
equipartition  of  the  turbulent  energy  remains  valid  for  [15],  [16],  and  [7].  It  appears  here  that 
stream-wise  fluctuations  contribute  the  most  to  the  average  turbulence  profile;  followed  by  the 
vertical  velocity  perturbations  and  the  lateral  velocity  perturbations.  As  the  distance  from  the  x- 
axis  increases,  stream-wise  fluctuations  decrease  quickly,  whereas  both  lateral  and  vertical  velocity 
perturbations  appear  not  to  decay  until  well  away  from  the  x-axis.  If  there  is  this  anisotropy  in 
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the  perturbation  profiles,  it  is  perhaps  natural  to  seek  experimental  data  which  is  not  azimuthally 
averaged. 

For  this  reason,  a  second  comparison  is  made  against  an  available  experimental  data  set  ( y  —  0) 
produced  by  Spedding  at  USC.  Because  of  the  experimental  techniques  used  to  collect  the  data, 
specific  planes  must  be  chosen  to  investigate  the  flowfield.  Thus,  the  velocity  field  in  the  vertical 
plane  is  available  for  comparison  (i.e.  Ux{z),  u'(z)  and  wf(z)).  The  perturbation  quantities  are 
defined  as: 

u'(z)  =  [u(x,t)-Ux(z)]x 
v'{z)  =  [t<(£,  t)]x 
w'{z)  =  [ui(f,f)]t 
Ux(z)  =  [u(x,  t)]x 

where  Ux(z)  is  now  the  mean  profile  in  the  vertical  slice,  the  subscript  x  indicates  spatial  averaging 
in  that  direction,  and  there  is  no  azimuthal  averaging  for  any  flow'  variable.  The  general  pattern  of 
comparison  from  before  is  continued  with  the  experimental  data.  Initial  comparisons  are  made  for 
the  predicted  mean  steam- wise  velocity  profile  against  the  experimented  profile  as  seen  in  Figure 
4.30.  Again,  without  the  azimuthal  averaging  the  DES  simulation  appears  to  predict  the  mean 
profile  with  some  degree  of  accuracy. 

Additional  comparisons  are  made  for  the  available  perturbation  velocities,  v!  and  wf.  The 
simulations  prediction  of  v!  against  experimental  data  is  seen  in  Figure  4.31.  The  model  appears 
to  over-predict  the  turbulent  fluctuations  in  the  ’’core”  (|z|  <  1)  of  the  wake  region.  However,  the 
model  still  captures  the  drop  in  fluctuations  near  the  center  of  the  core  at  z  =  0  for  v!  seen  in 
earlier  comparisons.  For  the  most  part,  the  stream-wise  fluctuations  are  located  within  the  core 
region.  Although  there  is  some  minor  diffusion  of  stream-wise  fluctuation  seen  out  of  the  core  for 
z  <  —  1,  perturbations  appear  fully  decayed  at  z  =  —  2. 

Remaining  comparisons  against  the  data  of  Spedding  may  then  be  made  for  the  vertical  velocity 
fluctuations,  it/.  In  the  Guassian  profile,  there  was  some  apparent  diffusion  out  of  the  turbulent 
core  for  w'  as  seen  in  4.29,  and  when  compared  to  experimental  data  in  the  vertical  plane  in  Figure 
4.32,  this  diffusion  of  w f  out  of  the  core  becomes  even  more  obvious.  The  SA-DES  model  appears 
to  capture  the  turbulent  levels  contained  within  the  core  \z\  <  —1  reasonably  well.  Once  again, 
the  reduction  of  turbulence  near  the  center  of  the  core  at  z  =  0  is  seen.  Major  levels  of  vertical 
perturbation  velocities  are  seen  outside  of  the  core  and  decay  of  these  perturbations  is  very  slow 
compared  to  the  stream-wise  fluctuation  profile  of  Figure  4.31.  This  region  of  high  vertical  velocity 
shall  be  referred  to  as  the  "tails”  of  the  vertical  perturbation  velocity  profiles. 

This  appears  to  indicate  some  issue  with  the  SA-DES  modeling  as  the  simulation  trend  does 
not  appear  to  correctly  follow  the  experimental  data  trend  produced  by  Spedding.  In  light  of 
Section  2.1.5,  which  was  in  turn  guided  by  conversations  with  other  investigators  in  this  field,  the 
suggestion  is  made  that  internal  waves  are,  at  this  point,  a  main  suspect  in  the  cause  of  these  w ' 
tails  due  to  the  reflection  of  the  internal  waves  off  of  the  domain  boundary. 

To  assess  this  suggestion,  the  method  of  Section  2.1.5,  a  zonal  DES  damper,  was  utilized  to 
investigate  the  possibility  of  a  quick  fix  to  this  behavior.  An  example  of  results  of  this  attempt 
are  seen  in  Figure  4.33.  It  is  clear  by  comparing  Figures  4.32  and  4.33  that  some  degree  of  w ' 
is  removed  from  outside  the  core  region.  However,  it  also  appears  that  now  fluctuations  have 
increased  inside  the  core  region.  This  result  gives  indication  that  some  boundary  treatment  may 
be  beneficial  to  the  accuracy  of  the  numerical  simulation.  Because  of  the  apparent  failure  of  the 
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SA-DES  RANS  damper,  a  more  recognized  method  is  pursued  in  the  form  of  the  artificial  viscous 
dampers  of  Section  2.1.5. 

Viscous  damping  to  the  w*  equation  is  investigated.  Viscous  dampers  are  applied  to  the  region 
within  the  computational  domain  where  no  data  collection  is  desired.  An  example  region  of  viscous 
damping  is  shown  in  Figure  4.34,  although  the  damping  region  may  be  applied  as  needed  but  guided 
by  the  discussion  of  Section  2.1.5.  An  example  effect  of  the  addition  of  the  viscous  damping  term 
is  shown  in  Figure  4.35.  Here,  the  perturbation  velocities  are  presented  in  the  azimuthal  average 
fashion  to  include  a  more  three  dimensional  view  of  the  effects  of  the  damper.  Clearly,  energy  has 
been  added  to  the  core  of  the  wake  while  the  ”  tails”  of  the  perturbation  velocity  have  not  even 
decayed  as  in  the  SA-DES  RANS  case.  This  trend  was  observed  for  a  wide  range  0(0.1)-0(100)  of 
the  i/(x)  coefficient  seen  in  Eq.  (2.16)  for  a  variety  of  damping  region  shapes. 

These  attempts  at  boundary  conditions  leave  the  question  of  the  origin  of  the  ’’tails”  unanswered 
as  of  yet.  Further  discussions  with  researchers  in  the  field  of  stratified  turbulence  have  suggested 
that  grid  resolution  is  perhaps  an  additional  reason  for  the  poor  comparison  with  experiments. 
Qualitative  inspection  compared  to  the  grid  resolution  may  be  compared  in  Figure  4.34.  The 
turbulent  regime  investigations  are  left  with  four  potential  ’’tail”  inducing  culprits:  the  artificial 
domain,  the  damper  type,  the  grid  resolution,  and  the  turbulence  model  itself. 
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Chapter  4 


Future  Work 


4.1  Laminar  Regime 

There  is  quite  a  possibility  that  the  impact  of  the  Schmidt  number  on  the  flow  around  the  sphere 
on  a  sphere  traveling  horizontally  through  a  stratified  fluid  has  not  yet  been  investigated  completely. 
This  problem  may  be  particularly  difficult  for  experimentalists  to  accomplish:  a  water  tank  is  not 
set  up  to  work  with  air,  and  a  wind-tunnel  is  not  particularly  adept  at  investigations  using  water 
as  its  working  fluid.  This  puts  the  current  investigation  in  a  unique  position  to  possibly  expand 
knowledge  in  this  niche  area. 

4.2  Turbulent  Regime 

A  two-tiered  approach  is  planned: 

First,  DNS  parallelized  computations  will  be  made  at  a  Reynolds  number  of  5000  across  a  range 
of  Froude  number.  As  shown  in  Section  3.2,  comparison  with  data  from  USC  is  available  as  a 
benchmark  for  computational  agreement.  This  approach  has  a  major  advantage  over  other  proposed 
approaches,  despite  its  high  computational  cost,  in  that  there  is  no  reliance  on  any  turbulence 
modeling.  Thus,  a  degree  of  complexity  is  removed  from  the  numerical  investigation,  which  will 
allow  the  first  tier  of  the  proposed  work  to  focus  on  the  underlying  physics  of  the  problem. 

A  second  benefit  of  the  DNS  approach  is  that  because  resolution  requirements  are  already  high, 
the  question  of  internal  wave  reflection  either  due  to  grid  coarsening  or  the  limitations  of  the  finite 
boundary  may  be  alleviated.  If  the  grid  is  found  to  be  sufficiently  fine,  but  reflections  are  still 
encountered  due  to  the  location  of  the  domain  boundary,  then  the  resolution  requirements  will  still 
aid  in  implementation  of  the  boundary  treatments  as  discussed  in  Section  2.1.5. 

The  second  step  of  the  two-tiered  approach  is  to  return  to  utilization  of  turbulence  modeling.  For 
all  its  benefits,  DNS  has  an  overwhelming  drawback:  with  current  world-wide  computational  re¬ 
sources,  DNS  is  not  a  tool  that  can  be  applied  for  geophysical  flows,  or  even  many  engineering-based 
flows.  Ultimately,  a  viable  turbulence  model  must  be  located  within  the  engineering  community  for 
the  simulation  of  the  turbulent  near-wake  of  flow  behind  a  sphere  traveling  horizontally  through  a 
stratified  fluid. 
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The  choice  of  turbulence  model  with  which  to  continue  investigation  is  currently  under  consid¬ 
eration.  The  Spalart  Allmaras  Detached  Eddy  Simulation  (SA-DES)  method  as  shown  in  Section 
3.2  has  not  produced  necessarily  successful  results.  However,  a  proclamation  of  the  SADES  method 
inappropriate  model. for  near-wake  flows  may  be  a  mistake.  Indeed,  its  ’’ease”  in  implementation, 
especially  for  pre-written  SA-RANS  based  codes,  makes  it  highly  tempting  as  a  turbulence  model. 
However,  investigations  utilizing  SA-DES  may  prove  as  difficult  as  investigations  with  any  other 
types  of  turbulence  model.  A  second  possible  choice  for  turbulence  modeling  is  Large  Eddy  Simu¬ 
lation  using  the  Truncated  Navier-Stokes  (LES-TNS)  method  (e.g.  Tantikul  &  Domaradzki  2010 
[53]).  This  has  an  advantage  over  the  SA-DES  method  in  that  USC  houses  the  resident  experts  in 
the  LES-TNS  method. 
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Figures 


Figure  4.1:  Sketch  of  Physical  Problem 


Figure  4.2:  Example  of  0-0  grid  topology 
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Figure  4.3:  Sketch  of  Body-Fitted  Coordinate  System 
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Figure  4.4:  Sketch  of  Computational  Grid 
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Figure  4.5:  Figure  2  of  Lin  et  al.  92  [32],  (Fi  =  Fd);  reproduced  here  for  convenience. 
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Figure  4.6:  Figure  3  of  Lin  et  al.  92  [32],  (Fi  =  FVi);  reproduced  here  for  convenience. 
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Figure  4.7:  Re=200,  Sc=l:  A  Co  vs.  1/  (2  Fd).  (□)  Orr  ’10,  (—  with  o)  Hanzaki  ’88,  (—  •  •—  with 
A)  Lofquist  &  Purtell  '84,  (—  with  V)  Mason  77. 
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Figure  4.8:  Re=200,  Sc=l:  Surface  pressure  coefficient,  Cp ,  vs.  angle,  6 ,  in  the  (-)  vertical(z  >  0) 
and  (—  ■  — )  horizontal  (y  <  0)  plane.  0  is  taken  from  the  negative  x-axis.  (a)  Fd=oc  (b)  Fd=l  (c) 
Fd=0.5  (d)  Fd=0.25  (e)  Fd=0.125 
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Figure  4.9:  Re=200,  Sc=l,  Fd=0.125:  Drag  coefficient  time  series.  —  Cdx)  —  •  —  —  Cdz. 


Figure  4.10:  Re=200,  Sc=l:  Drag  coefficient  power  spectra.  -  Cdx,  —  •  —  Cdv, - Cdz- 


31 
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Figure  4.11:  Re=200,  Sc=l,  Fd=0.125:  Vertical  vorticity  contours  in  the  horizontal  plane  at  z=0. 
Contours  of  ujz  are  leveled  between  ±1  in  increments  of  0.1 
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Figure  4.12:  Re=200,  Sc=l:  Surface  skin  friction  coefficient,  Cj,  vs.  angle,  6,  in  the  (— )  vertical(z  > 
0)  and  (—  •  — )  horizontal  ( y  <  0)  plane.  6  is  taken  from  the  negative  x-axis.  (a)  Fd=oo  (b)  Fd=l 
(c)  Fd=0.5  (d)  Fd=0.25  (e)  Fd=0.125 
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Figure  4.13:  Re=200  Separation  point  location  vs.  1/  (2  Fd).  Orr  ’10  (Sc=l):  (v)  vertical,  (0) 
horizontal.  Hanazaki  ’88  (Sc=oc):  (— )  vertical,  ( — )  horizontal. 
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Figure  4.14:  Re=200,  Sc=l:  Streamlines  on  the  lee  side  of  the  sphere,  (a)  Fd=oc,  (b)  Fd=1.0,(c) 
Fd=0.5,  (d)  Fd=n.25,  (e)  Fd=0.125 
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Figure  4.15:  Re=200,  Sc=l:  Isopycnal  lines.  Levels  of  density  axe  chosen  between  ±1.5  in  incre¬ 
ments  of  0.075.  (a)  Fd=l  (b)  Fd=0.50  (c)  Fd=0.25  (d)  Fd=0.125 
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(c)  (d) 


Figure  4.16:  Re=200,  Sc=l:  Contours  of  vertical  velocity,  w.  Levels  are  between  ±.01  in  increments 
of  0.001.  (a)  Fd=l,  (b)  Fd=0.50,  (c)  Fd=0  25,  (d)  Fd=0.125 
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Figure  4.17:  Re=1000,  Sc=l:  A  Co  vs.  1/  (2  Fd).  (□)  Orr  ’10,  (-  with  o)  Hanzaki  ’88,  (-  • 
with  A)  Lofquist  &  Purtell  ’84,  (-  with  V)  Mason  77. 
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Figure  4.18:  Rc=1000,  Sc=l:  Time-averaged  surface  coefficient  of  pressure,  Cp,  vs.  angle,  9 ,  in 
the  (  — )  verticals  >  0)  and  (—  •  — )  horizontal  ( y  <  0)  plane.  9  is  taken  from  the  negative  a:-axis. 
(a)  Fd=oc  (b)  Fd=2  (c)  Fd=0.50  (e)  Fd=0.25 
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Figure  4.19:  Re-1000,  Sc=l:  Surface  skin  friction  coefficient,  C/,  vs.  angle,  ft,  in  the  (-) 
vertical(z  >  0)  and  (—  •  — )  horizontal  ( y  <  0)  plane,  ft  is  taken  from  the  negative  x-axis.  (a) 
Fd=oo  (b)  Fd=2  (c)  Fd=0.50  (e)  Fd=0.25 
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Figure  4.20:  Re=1000,  Sc=l:  Separation  point  location  vs.  1/  (2  Fd).  (v)  vertical,  (<C>)  horizontal. 
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(a) 


(b) 


(c)  (d) 

Figure  4.21:  Re=1000,  Sc=l:  Instantaneous  vectors  in  vertical  plane  ( z  >  0)  on  lee  side  of  sphere, 
(a)  Fd=oo  (b)  Fd=2  (c)  Fd=0.50  (d)  Fd=0.25 
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Figure  4.22:  Re=1000,  Sc=l,  Fd=0.25:  Drag  coefficient  time  series.  —  Cdv,  —  •  —  Cdy,  —  Cdz- 
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Figure  4.23:  Re=1000,  Sc=l,  Fd=0.25:  Drag  coefficient  power  spectra.  —  CpT,  -  *  —  Cdu, 
cDt- 


(a)  (b) 


Figure  4.24:  Re=1000,  Sc=l,  Fd=2:  Time  history  of  data  signal  at  x/Dc^2.  (a)  w  (b)  pf 


(a)  (b) 

Figure  4.25:  Re=1000,  Sc=l,  Fd=2:  Power  spectra  of  data  signal  at  x/D^2.  (a)  w  (b)  p' 


43 
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(c) 


Figure  4.26:  Re=1000,  Sc=l:  Contours  of  vertical  velocity,  w.  Levels  are  between  ±.01  in  incre¬ 
ments  of  0.001.  (a)  Fd=2,  (b)  Fd=0.50,  (c)  Fd=0.25 
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Figure  4.27:  Re=5000,  Sc=l,  Fd=2:  Azimuthally  averaged  stream-wise  velocity  profile.  —  *  -  Orr 
’10,  -  Diamessis  ’OS.  (i.c.) 


r/D 


Figure  4.28:  Re=5000,  Sc=l,  Fd=2:  Azimuthally  and  spatially  averaged  (uj)  velocity  profile  at 
(1  /Fd)(x/D)  =3.  ( - )  Orr  ’10,  (-)  Diamessis  ’05.  (i.c.) 
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Figure  4.29:  Re=5000,  Sc=l,  Fd=2:  Azimuthally  and  spatially  averaged  u\  velocity  profile  at 
(1  /Fd)(x/D)  =  3.  ({u'),  {?/),  (w'))  Orr  '10,  (  — )  Diamessis  '05.  (i.c.) 


Slice  along  y  =  0 


Figure  4.30:  Re=5000,  Sc=l,  Fd=2:  Spatially  averaged  ux  velocity  profile  in  vertical  center-plane 
(y  =  0)  at  (1  /Fd)(x/D)  =  3.  (-  •  -)  Orr  ’10,  (•)  Spedding 
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Slie*  along  y  =  0 


Figure  4.31:  Re=5000,  Sc=l,  Fd=2:  Spatially  averaged  uf  velocity  profile  in  vertical  center-plane 
( y  =  0)  at  (1  /Fcl)(x/D)  =3.  (-  •  -)  Orr  '10,  (•)  Spedding 


Slice  along  y  =  0 


Figure  4.32:  Re=5000,  Sc=l,  Fd=2:  Spatially  averaged  wf  velocity  profile  in  vertical  center-plane 
( y  =  0)  at  (1  /Fd)(x/D)  =  3.  (—  •  — )  Orr  ’10,  (•)  Spedding 
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Figure  4.33:  Re=5000,  Sc=l,  Fd=2:  Example  of  SA-DES  RANS  damper  affect  on  w'  velocity 
profile  in  vertical  center-plane  (y  =  0)  at  (1  /Fd)(x/D)  =  3.  (—  •  — )  Orr  ?10,  (•)  Spedding 


Figure  4.34:  Example  of  viscous  damping  layer  region  in  the  vertical  plane  (z  >  0);  0  indicates  no 
damping,  1  indeates  full  damping.  Mesh  density  visible  in  vertical  plane  (z  <  0)  for  qualitative 
comparison. 
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Figure  4.35:  Re=5000,  Sc=l,  Fd=2:  Example  of  viscous  damper  affect  on  vJ  velocity  profile  in 
vertical  center-plane  (y  =  0)  at  (1  /Fd)(x/D)  =  3.  (—  •  — )  Orr  ’10,  (•)  Spedding 


19 


Bibliography 


[1]  A.  Arnone,  M.S,  Liou,  and  L.A.  Povinelli.  Integration  of  Navier-Stokes  equations  using  dual 
time  stepping  and  a  multigrid  method.  AIAA  Journal ,  33:985-990,  1995. 

[2]  P.G.  Baines.  Upstream  blocking  and  airflow  over  mountains.  Annual  Review  of  Fluid  Mechan¬ 
ics, ,  19:75-97,  1987. 

[3]  P.M  Bevilaqua  and  P.S.  Lykoudis.  Tubulence  memory  in  self- preserving  wakes.  Journal  of 
Fluid  Mechanics ,  89:589-  606,  1978. 

[4]  P.  Bonneton,  J.M.  Chomaz,  E.J.  Hopfinger,  and  M.  Perrier.  The  structure  of  the  turbulent 
wake  and  the  random  internal  wave  field  generated  by  a  moving  sphere  in  a  stratified  fluid. 
Dynamics  of  Atmospheres  and  Oceans ,  23:299-308,  1996. 

[5]  M.  Bonnier,  P.  Bonneton,  and  O.  Eiff.  Far-wake  of  a  sphere  in  a  stably  stratified  fluid: 
Characterization  of  the  vortex  structures.  Applied  Scientific  Research ,  59:269  281,  1998. 

[6]  M.  Bonnier,  O.  Eiff,  and  P.  Bonneton.  On  the  density  structure  of  far-wake  vortices  in  a 
stratified  fluid.  Dynamics  of  Atmospheres  and  Oceans ,  31:117-137,  2000. 

[7]  K.A.  Brucker  and  S.  Sarkar.  A  comparative  study  of  self-propelled  and  towed  wakes  in  a 
stratified  fluid.  Journal  of  Fluid  Mechanics ,  652:373-404,  2010. 

[8]  J.  D.  Chashechkin  and  H.  J.  Sysoeva.  Fine  structure  and  symmetry  of  the  wake  past  a  sphere 
in  a  stratified  liquid.  Proceedings  of  17th  Session  of  D.S.H.C.,  1:10-1,  1988. 

[9]  J.M.  Chomaz,  P.  Bonneton,  and  E.J.  Hopfinger.  The  structure  of  the  near  wake  of  a  sphere 
moving  horizontally  in  a  stratified  fluid.  Journal  of  Fluid  Mechanics ,  254:1-21,  1993. 

[10]  J.M.  Chomaz,  P.  Bonneton,  M.  Perrier,  and  E.J.  Hopfinger.  Froude  number  dependence  of 
the  flow  separation  line  on  a  sphere  towed  in  a  stratified  fluid.  Physics  of  Fluids,  4:254  258, 
1992. 

[11]  G.S.  Constant inescu,  M.  Chapelet,  and  K.D.  Squires.  Turbulence  modeling  applied  to  flow 
over  a  sphere.  AIAA  Journal ,  41:1733-1742,  2003. 

[12]  G.S.  Constantineseu  and  V.C.  Patel.  Numerical  model  for  simulation  of  pump-intake  flow  and 
vortices.  Journal  of  Hydraulic  Engineering ,  124:123  134,  1998. 

[13]  G.S.  Constantineseu  and  K.D.  Squires.  LES  and  DES  investigations  of  turbulent  flow  over  a 
sphere  at  re  =  10,000.  Flow ,  Turbulence  and  Combustion ,  70:267-298,  2003. 


50 


[14]  G.S.  Constantinescu  and  K.D.  Squires.  Numerical  investigations  of  flow  over  a  sphere  in  the 
subcritical  and  supercritical  regimes.  Physics  of  Fluids,  16:1449-1466,  2004. 

[15]  P.J.  Diamessis,  J. A.  Domaradzki,  and  J.S.  Hesthaven.  A  spectral  multidomain  penalty  method 
model  for  the  simulation  of  high  reynolds  number  localized  incompressible  stratified  turbulence. 
Journal  of  Computational  Physics ,  202:298-322,  2005. 

[16]  D.G.  Dommermuth,  Rottman  J.W.,  G.E.  Innis,  and  E.A.  Novikov.  Numerical  simulation  of 
the  wake  of  a  towed  sphere  in  weakly  stratified  fluid.  Journal  of  Fluid  Mechanics ,  473:83-101, 
2002. 

[17]  P.G.  Drazin.  On  the  steady  flow  of  a  fluid  of  variable  density  over  an  obstacle.  Tellus ,  8:239- 
251,  1961. 

[18]  D.  Givoli.  High-order  local  non- reflecting  boundary  conditions:  a  review.  Wave  Motion , 
39:319-326,  2004. 

[19]  M.P.  Gresho  and  R.L.  Sani.  On  pressure  boundary  conditions  for  the  incompressible  Navier- 
Stokes  equations.  International  Journal  for  Numerical  Methods  in  Fluids ,  7:1111-1145,  1987. 

[20]  G.  Haller.  Exact  theory  of  unsteady  separation  for  two-dimensional  flows.  Journal  of  Fluid 
Mechanics ,  512:257-311,  2004. 

[21]  H.  Hanazaki.  A  numerical  study  of  three-dimensional  stratified  flow  past  a  sphere.  Journal  of 
Fluid  Mechanics ,  192:393-419,  1988. 

[22]  H  Hanazaki,  K.  Konishi,  and  T.  Okamura.  Schmid t-number  effects  on  the  flow  past  a  sphere 
moving  vertically  in  a  stratified  diffusive  fluid.  Physics  of  Fluids ,  21:1-8,  2009. 

[23]  E.J.  Hopfinger,  J.-B.  Flor,  J.-M.  Comaz,  and  P.  Bonneton.  Internal  waves  generated  bv  a 
moving  sphere  and  its  wake  in  a  stratified  fluid.  Experiments  in  Fluids ,  11:255-261,  1991. 

[24]  M.  Israeli  and  S.  A.  Orszag.  Approximation  of  radiation  boundary  conditions.  Journal  of 
Computational  Physics ,  41:115  135,  1981. 

[25]  T.A.  Johnson  and  V.C.  Patel.  Flow  past  a  sphere  up  to  a  Reynolds  number  of  300.  PhD  thesis, 
University  of  Iowa,  1999. 

[26]  Pijush  Kundu  and  Ira  Cohen.  Fluid  Mechanics ,  Third  Edition.  Elsevier  Inc.,  2004. 

[27]  S.  Lee.  A  numerical  study  of  the  unsteady  wake  behind  a  sphere  in  a  uniform  flow  at  moderate 
reynolds  numbers.  Computers  &  Fluids,  pages  639  667,  2000. 

[28]  James  Light  hill .  Waves  in  Fluids.  Cambridge  University  Press,  New  York,  NY,  1978. 

[29]  D.K.  Lilly.  The  representation  of  small-scale  turbulence  in  numerical  simulation  experiments. 
In  Proc.  IBM  Sci.  Comput.  Symp.  Environ.  Sci. ,  pages  195-210,  White  Plains,  N.Y.,  1967. 
IBM  Data  Process.  Div. 

[30]  J.T.  Lin  and  Y.H.  Pan.  Wakes  in  stratified  fluids.  Annual  Review  of  Fluid  Mechanics ,  1 1:317— 
338,  1979. 


51 


[31]  Q.  Lin,  D.L.  Boyer,  and  H.J.S.  Fernando.  Turbulent  wakes  of  linearly  stratified  flow  past  a 
sphere.  Physics  of  Fluids  A ,  4:1687-1696,  1992. 

[32]  Q.  Lin,  W.R.  Lindberg,  D.L.  Boyer,  and  H.J.S.  Fernando.  Stratified  flow  past  a  sphere.  Journal 
of  Fluid  Mechanics,  240:315-  354,  1992. 

[33]  Vladimir  Liseikin.  Grid  Generation  Methods.  Springer,  1999. 

[34]  K.E.  Lofquist  and  L.P.  Purtell.  Drag  on  a  sphere  moving  horizontally  through  a  stratified 
liquid.  Journal  of  Fluid  Mechanics,  148:271-284,  1984. 

[35]  P.  Meunier  and  G.R.  Spedding.  A  loss  of  memory  in  stratified  momentum  wak^s.  Physics  of 
Fluids ,  16(2):298  305,  February  2004. 

[36]  P.  Meunier  and  G.R.  Spedding.  Stratified  propelled  wakes.  Journal  of  Fluid  Mechanics, 
552:229-256,  2006. 

[37]  N.V.  Nikitin,  F.  Nicoud,  B.  Wasistho,  K.D.  Squires,  and  P.R.  Spalart.  An  approach  to  wall 
modeling  in  large-eddy  simulations.  Physics  of  Fluids,  12(7):162£  1632,  July  2000. 

[38]  H.-P.  Pao  and  T.W.  Kao.  Vortex  structure  in  the  wake  of  a  sphere.  The  Physics  of  Fluids, 
20:187-191,  1977. 

[39]  L.  Prandtl.  Uber  fliissigkcitsbewegung  bei  sehr  kleiner  reibung.  Int.  Math.  Kongr.  Heidelberg , 
pages  484  451,  1904. 

[40]  H.  Schlichting.  Boundary- Layer  Theory.  McGraw-Hill  Book  Company,  7  edition,  1979. 

[41]  A.  H.  Schooley  and  R.  W.  Stewart.  Experiments  with  a  self-propelled  body  submerged  in  a 
fluid  with  a  vertical  density  gradient.  Journal  of  Fluid  Mechanics,  9:83-96,  1963. 

42]  W.R.  Sears  and  D.P.  Telionis.  Boundary-layer  separation  in  unsteady  flow.  SIAM  Journal  on 
Applied  Mathematics ,  28:215-235,  1975. 

[43]  P.A.  Sheppard.  Airflow  over  mountains.  Quarterly  Journal  of  the  Royal  Meteorological  Society , 
82:528-529,  1956. 

[44]  M.  Shur,  P.R.  Spalart,  M.  Strelets,  and  A.  Travin.  Detached-eddv  simulation  of  an  airfoil  at 
high  angle  of  attack.  In  Proceedings  of  4th  international  symposium  on  engineering  turbulence 
modelling  and  measurements.  Elsevier,  May  1999. 

[45]  F.  Sortiropoulos.  A  primitive  variable  method  for  the  solution  of  external  and  internal  incom¬ 
pressible  flow-fields.  PhD  thesis,  University  of  Cincinnati,  1991. 

[46]  F.  Sortiropoulos  and  S.  Abdallah.  The  discrete  continuity  equation  in  primitive  variable  solu¬ 
tions  of  incompressible  flow.  Journal  of  Computational  Physics,  95:212  227,  1991. 

[47]  P.R.  Spalart  and  S.R.  Allmaras.  A  one-equation  model  for  aerodynamic  flows.  La  Recherch 
Aerospatiale,  1:1,  1994. 


52 


[48]  P.R.  Spalart,  W-H.  Jon,  M.  Strelets,  and  S.R.  Allmaras.  Comments  on  the  feasibility  of  les 
for  wings,  and  on  a  hybrid  rans/les  approach.  In  1st  AFOSR  Int.  Conf.  on  DNS/LES.  1st 
AFOSR  Int.  Conf.  on  DNS/LES,  Aug.  4-8  1997. 

[49]  G.R.  Spedding.  The  evolution  of  initially  turbulent  bluff-body  wakes  at  high  internal  froude 
number.  Journal  of  Fluid  Mechanics ,  337:283*301,  1997. 

[50]  M.  Strelets.  Detached  eddy  simulation  of  massively  separated  flows.  In  39th  AIAA  Aerospace 
Sciences  Meeting  and  Exhibit  39th  AIAA  Aerospace  Sciences  Meeting  and  Exhibit,  Jan.  8-11 
2001. 

[51]  E.  Ya.  Sysoeva  and  Yu.  D.  Chashechkin.  Vortex  structure  of  a  wake  behind  a  sphere  in  a 
stratified  fluid.  J.  Appl.  Mech.  Theor.  Phys.  ( Novosibirsk ),  2:190-196,  1986. 

[52]  S.  Taneda.  Experimental  investigation  of  the  wake  behind  a  sphere  at  low  reynolds  numbers. 
Journal  of  the  Physical  Society  of  Japan,  11:1104-1108,  1956. 

[53]  T.  Tantikul  and  J.A.  Domaradzki.  Large  eddy  simulations  using  truncated  Navier- Stokes 
equations  with  the  automatic  filtering  criterion.  Journal  of  Turbulence ,  11,  2010. 

[54]  L.L.  Van  Dommelen  and  S.J.  Cowley.  On  the  lagrangian  description  of  unsteady  boundary- 
layer  separation,  parti,  general  theory.  Journal  of  Fluid  Mechanics ,  210:593-626,  1990. 

[55]  D.  Wilcox.  Turbulence  Modeling  for  CFD.  DCW  Industries,  Inc.,  1998. 


53 


Downloaded  By:  [USC  University  of  Southern  California]  At:  21:56  28  July  2010 


Journal  of  Turbulence 
Vol.  11,  No.  21,  2010,  1-24 


©Taylor  &  Francis 

teyfcwfc  Franc*  Oft* 


Large  eddy  simulations  using  truncated  Navier-Stokes  equations  with 
the  automatic  filtering  criterion 

T.  Tantikul*  and  J.A.  Domaradzki 

Department  of  Aerospace  and  Mechanical  Engineering,  University  of  Southern  California, 

Los  Angeles,  CA,  90089-1191,  USA 

{Received  28  November  2009;  final  version  received  30  April  2010) 

We  propose  a  large  eddy  simulation  (LES)  technique  based  on  the  previously  developed 
Trancated  Navier-Stokes  (TNS)  method.  In  TNS,  the  Navier-Stokes  equations  are 
solved  through  a  sequence  of  direct  numerical  simulation  runs  and  a  penodic  processing 
of  small  scales  to  provide  the  necessary  dissipation.  In  the  simplest  case,  the  processing 
is  accomplished  by  filtering  the  turbulent  fields  with  a  properly  chosen  filter  In  the 
previous  work,  the  period  for  processing  was  selected  in  advance  for  each  case  using 
heuristic  arguments  validated  by  trial  and  error  In  this  work,  we  develop  a  criterion  that 
automates  the  selection  of  a  time  instant  in  simulations  when  the  processing  occurs. 

The  criterion  is  based  on  the  relationship  between  the  energy  of  the  flow  field  and  the 
energy  of  the  same  field  filtered  with  the  chosen  filter.  The  procedure  is  tested  in  LES 
of  the  turbulent  channel  flow  performed  at  various  Reynolds  numbers  and  in  domains 
of  different  sizes  for  which  Direct  Numerical  Simulations  (DNS)  and  experimental  data 
are  available  for  comparisons. 

Keywords:  large  eddy  simulations;  truncated  Navier-Stokes  equations;  explicit  filter¬ 
ing;  automatic  filtering  criterion 


1.  Introduction 

Several  classifications  of  the  subgrid-scale  (SGS)  models  for  large  eddy  simulations  (LES) 
have  been  proposed  in  the  literature  on  the  subject.  In  a  rev  iew  paper  of  Domaradzki  and 
Adams  [1],  the  SGS  models  are  divided  into  two  general  categories:  the  traditional  models 
that  use  the  explicit  expressions  for  the  SGS  terms,  the  SGS  stress  tensor  in  particular  and 
the  models  that  construct  the  unknown  primitive  variables,  such  as  velocity,  in  order  to 
use  these  variables  to  compute  the  SGS  terms  directly  from  the  definitions.  The  traditional 
SGS  modeling  approaches  can  be  subdivided  into  three  groups:  the  eddy  viscosity  models, 
the  similarity  models  and  the  mixed  models.  The  examples  of  the  models  in  the  other 
group  are  the  velocity  estimation  model  proposed  by  Domaradzki  et  al.  [2,3]  and  Stolz 
et  al/s  [4]  approximate  deconvolution  model  (ADM).  Kosovic  et  al.  [5],  characterized  by 
the  use  of  explicit  filtering,  divided  SGS  models  into  three  groups.  First  are  the  models 
that  use  an  explicit  SGS  expression  but  do  not  employ  explicit  filtering;  for  instance,  the 
classical  Smagorinsky  eddy  viscosity  model.  Second  are  the  models  that  use  an  explicit  SGS 
expression  and  employ  explicit  filtering,  c.g.  the  dynamic  Smagorinsky  model.  The  third 
category  are  the  so-called  implicit  models  because  equations  of  motion  are  solved  without 
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neither  explicit  filtering  nor  explicit  SGS  terms.  Such  approaches  rely  on  the  properties  of  the 
numerical  scheme  to  provide  implicit  dissipation  and  are  known  as  the  Monotone  Integrated 
Large  Eddy  Simulations  (MILES)  or  the  Implicit  Large  Eddy  Simulations  (ILES).  An  in 
depth  review  of  the  ILES  is  given  in  the  monograph  edited  by  Grinstein  et  al.  [6].  For  high- 
ordcr,  nondissipativc  numcncal  schemes,  this  classification  can  be  extended  to  include 
models  that  do  not  use  the  explicit  SGS  terms  but  employ  the  explicit  filtering  to  model  the 
effects  of  the  SGS  dissipation.  The  modeling  approach  discussed  in  this  paper  belongs  to  this 
group  and  originated  from  the  SGS  velocity  estimation  model  proposed  by  Domaradzki  and 
Saiki  [2].  The  model  was  developed  in  spectral  space  representation  and  later  implemented 
in  the  physical  space  by  Domaradzki  and  Loh  [7]. 

The  original  velocity  estimation  model  encountered  difficulties  in  the  LES  of  high 
Reynolds  number  flows  because  it  did  not  produce  enough  SGS  dissipation.  To  address 
this  shortcoming,  Domaradzki  ct  al.  [3]  have  introduced  the  modification  of  the  velocity 
estimation  model  that  solves  LES  equations  in  parallel  with  the  Truncated  Navier-Stokes 
(TNS)  equations  on  a  finer  mesh;  TNS  provides  the  estimated  velocity  that  is  used  to 
compute  the  SGS  stress  tensor  needed  in  LES.  More  specifically,  the  TNS  equations  are 
just  Navier-Stokes  equations  discretized  on  a  mesh  by  a  factor  of  two  smaller  than  the  LES 
equations  of  interest.  The  estimated  velocity  to  be  advanced  in  TNS  is  obtained  in  two  steps. 
First,  the  velocity  from  LES  is  ‘de-filtered’,  providing  large-scale  velocity  represented  on 
a  coarse  LES  mesh.  Second,  the  small-scale  perturbation  velocity  is  created  by  nonlinear 
interaction  of  the  reconstructed  field.  The  total  estimated  velocity  is  a  sum  of  both  parts 
and  is  represented  on  a  finer  TNS  mesh.  The  estimated  velocity  is  advanced  in  time  using 
TNS  equations,  i.e.  using  a  Navier-Stokes  solver  on  the  fine  mesh.  The  velocity  field  from 
TNS  is  then  used  at  each  time-step  to  the  SGS  stress  model  directly  from  the  TNS  field,  and 
this  stress  model  is  used  to  advance  in  time  the  LES  equations.  This  process  of  two  parallel 
runs  is  continued  for  some  fixed  time  period  T  until  the  accumulation  of  energy  in  the  small 
scales  of  the  TNS  begins  to  affect  the  large  scales.  To  assure  the  accurate  large-scale  results, 
the  whole  process  has  to  be  reinitialized  every  time  period  T,  with  reduced  energy  in  the 
small-scale  region  obtained  by  a  low-pass  filter.  The  results  from  LES  for  homogeneous 
turbulence  and  turbulent  channel  flow  test  cases  are  in  a  very  good  agreement  with  the 
Direct  Numerical  Simulations  (DNS)  and  the  experimental  data.  The  disadvantage  of  this 
approach  is  that  the  procedure  is  quite  complicated  and  requires  more  computational  time 
compared  to  other  the  LES  modeling  approaches  because  additional  TNS  equations  must  be 
solved.  However,  a  significant  speed  up  of  computations  is  possible  for  a  simplified  version 
of  the  model  that  employs  only  the  TNS  equations  as  shown  by  Domaradzki  et  al.  [3]. 
Such  a  simplification  is  based  on  the  observation  that  the  TNS  velocity  on  the  fine  mesh 
already  contains  information  about  large  scales  of  the  flow  and  thus  the  LES  equations 
are  redundant.  Since  the  small  scales  of  TNS  fields  are  subjected  to  periodic  filtering,  the 
physical  meaning  is  ascribed  only  to  the  large  TNS  scales  given  on  the  coarse  LES  mesh. 
In  this  approach  the  estimated  small-scale  velocity  field  provides  the  SGS  effects  on  the 
large  scales  via  nonlinear  interactions  in  the  Navier-Stokes  equations.  The  approach  is 
less  complicated  in  terms  of  the  numerical  implementation  and  retains  the  same  quality  of 
the  results  as  the  original  method.  However,  the  user  still  needs  to  specify  the  time  period 
T  for  filtering  and  re-initialization  of  the  TNS  in  order  to  prevent  the  accumulation  of 
energy  in  the  small  scales  and  the  resulting  error  propagating  toward  the  large  scales  and 
contaminating  them.  Our  goal  in  this  work  is  to  investigate  a  procedure  that  determines  the 
time  period  T  automatically,  without  a  need  to  provide  its  value  as  a  parameter  by  the  user. 
Once  the  procedure  is  defined  in  the  next  section  it  will  be  subsequently  tested  in  LES  of 
the  turbulent  channel  flow  problem  in  the  following  sections. 
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The  explicit  filtering  as  a  LES  tool  has  been  employed  in  different  forms  previously, 
usually  as  a  method  to  control  numerical  instabilities.  Visbal  and  Rizzetta  [8],  [9]  studied 
the  application  of  a  high-order  finite  difference  method  to  the  LES  of  compressible  flows. 
They  used  fourth-  and  sixth-order  compact  differencing  schemes  for  spatial  discretization 
and  different  explicit  and  implicit  time  marching  methods.  The  simulations  were  performed 
with  the  unfiltered  Navicr-Stokes  equations  and  the  high-order,  Pade-type,  low-pass  spatial 
filter  was  employed  at  each  time-step  to  remove  spurious  high-frequency  modes  which 
arise  because  of  the  lack  of  the  dissipation  in  the  spatial  discretization.  It  was  found  that  the 
compact/filtering  scheme  performed  better  or  comparably  to  the  constant  coefficient  and 
dynamic  Smagorinsky  models  for  a  number  of  different  turbulent  flows.  Bogey  and  Bailly 
[10]  investigated  jet  flows  using  LES  with  explicit  filtering.  The  numerical  simulations 
were  performed  using  a  low-dissipation  numerical  scheme  with  the  explicit  fourth-order 
1 3-point  centered  finite  differences  for  spatial  discretizations  and  the  second-order  six-stage 
low-storage  Runge-Kutta  algorithm  for  time  integration.  The  explicit  selective/high-order 
filtering  was  employed  to  emulate  the  dissipative  effects  of  the  neglected  SGS  Interestingly, 
the  filtering  was  not  applied  at  every  time-step  but  every  second  or  third  time-step.  The 
results  from  the  simulations  with  the  filtering  applied  at  every  two  or  three  iterations  exhibit 
no  relationship  between  the  frequency  of  application  of  the  filtering  and  the  quality  of  the 
results.  When  spectral  methods,  which  have  negligible  numerical  dissipation,  are  employed 
to  simulate  high  Reynolds  number  flows,  explicit  filters  and  penalty  techniques  must  be 
used  to  stabilize  numerics,  e.g.  in  Diamessis  et  al.  [1 1],  and  in  Minguez  and  Pasquetti  [12]. 
Such  methods  are  known  as  the  stabilized  LES  and  use  filtering  at  each  time  step.  The 
numerical  dissipation  associated  with  one  of  these  methods  has  been  investigated  recently 
by  Diamessis  et  al.  [13].  In  using  such  techniques  one  must  recognize  that  achieving 
numerical  stability  does  not  necessarily  imply  physically  correct  results.  Therefore,  such 
techniques  must  always  be  validated  by  comparing  results  they  produce  with  experiments 
and  DNS  results. 

In  the  work  described  here,  we  use  spectral  methods  stabilized  by  the  weak,  explicit 
filter  applied  at  each  time  step  and  much  less  frequent  periodic  filtering  that  is  used  not  as  a 
stabilization  tool  but  primarily  as  an  SGS  modeling  tool  based  on  physical  considerations 
of  the  energy  transfer  and  the  spectra  in  developed  turbulence. 


2.  Filtering  and  filtering  criterion 

In  LES,  the  filtering  operation  is  applied  to  a  full  and  turbulent  velocity  field,  to  separate  it 
into  the  large,  resolved  scales  and  the  small  SGS.  Ideally,  the  filter  should  retain  complete 
information  about  the  large  scales  and  substantially  remove  or  attenuate  the  small  scales, 
allowing  the  user  to  focus  the  computational  power  on  the  large  scales.  Since  the  compu¬ 
tational  space  in  the  simulations  is  discrete,  filters  are  given  in  the  discrete  form  as  well. 
A  convenient  filter  with  desirable  properties  has  been  proposed  by  Stolz  et  al.  [4]  for  the 
Approximate  Deconvolution  Model  (ADM).  If  G  represents  the  primary  filter  used,  the 
filtering  operation  can  be  written  in  the  physical  space  as 

/OG 

G(x  —  x)u{x)dx  (1) 

-OG 

and  in  the  spectral  space  as 


u(k)  =  G(k)u(k ), 


(2) 
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where  u  is  the  one-dimensional  velocity  component,  u  is  the  filtered  velocity,  G(x)  is  the 
filter  kernel  in  the  physical  space  and  G(k)  is  the  filter  kernel  or  the  filter  transfer  function 
in  Fourier  space.  The  approximate  deconvolution  of  the  primary  filter  is  giv  en  as 

N 

Qh*G-'  =  £(/-G)v.  (3) 

u=0 


The  secondary  filter,  which  affects  only  the  small  resolved  scales,  is  constructed  in  ADM  as 
a  product  of  Qh  and  G  [4]  and  the  corresponding  filtered  velocity  is  expressed  as  follows: 

u  =  (QnG)*u.  (4) 

The  numerical  evaluation  of  the  filter  operation.  Equation  (1),  requires  the  quadrature 
rule.  We  select  the  box  filter  with  filter  width  A  as  the  primary  filter  and  the  trapezoidal 
integration  rule  and  A  =  2 h  for  physical  LES  scales.  For  this  choice,  the  one-dimensional 
transfer  function  is 


G(Jfc)  =  A(1  +  cos  kh).  (5) 

The  discrete  filter  on  equidistant  mesh  in  the  x  direction  is 

l 

«(•*/)  =  dgu(Xi)  +  ^2  d“u(Xi+i),  (6) 

/=— i,/#o 

where  d“  are  filter  coefficients  computed  by  Loh  et  al.  [14]  using  the  numerical  integration 
with  trapezoidal  rule  over  interval  A  spanning  three  neighboring  points  on  a  uniform  mesh, 
i.e.  A  =  2h,  where  A  is  the  filter  width.  The  filtering  in  the  simulations  is  done  sequentially 
in  each  Cartesian  direction.  For  the  nonuniform  mesh  in  the  vertical  direction,  the  discrete 
filter  has  the  same  form  as  given  in  Equation  (6)  but  the  filter  coefficients  are  different,  d”, 
with  the  explicit  formulas  given  in  Loh  et  al.  [14]. 

Jcanmart  and  Winckclmans  [15]  proposed  a  significant  simplification  of  the  filtering 
procedure  in  Equation  (4),  which  dispenses  entirely  with  the  the  Van  Cittert  approximate 
deconvolution  in  Equation  (3).  Indeed,  using  the  formula  for  the  sum  of  a  geometric 
sequence 


£ 


u=0 


1  —  rN+l 

— : - .  r  £  1, 

1  —  r 


(7) 


Equation  13)  is  expressed  formally  as 


Qn  = 


/  -  (/  -  C)  v+1 
G 


and  the  transfer  function  of  the  ADM  filter  in  Equation  (4)  becomes 


(8) 


QnG  =  /-  (/-  G)N+l, 


(9) 
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i.e.  the  ADM  filter  QsG  on  the  LHS  is  easily  implemented  through  the  recursive  application 
of  the  (/  —  G)  filter  on  the  RHS.  Also,  the  Equation  (5),  which  is  the  transfer  function  of 
the  stencil-3,  second-order  filter,  can  be  rearranged  as  follows: 


(10) 


which  immediately  leads  to  the  expression  for  the  transfer  function  for  the  ADM  filter  (9) 


(11) 


Equation  (1  1)  is  the  transfer  function  of  the  high-order  filter  obtained  from  the  recursive 
application  of  the  stencil-3,  low-order  filter  (/  —  G).  This  filter  can  be  applied  sequentially 
in  each  direction  in  three-dimensional  space.  In  particular,  the  1-D  filter  selected  in  this 
work  is  Q$G,  with  the  transfer  function  (1  -  sinn(kh/2)). 

In  a  simulation  with  general  finite  difference  schemes,  the  inability  to  represent  the 
high  wave  number  modes  accurately  results  in  the  undesired  dispersion  error.  This  error 
causes  unphysical  and  rapid  oscillations  in  the  marginally  resolved  regions,  which  may  lead 
to  numerical  instabilities  and  a  break  down  of  the  simulation.  To  eliminate  such  spurious 
modes  an  artificial  dissipation  through  additional  damping  terms  is  usually  employed.  The 
spectral  methods  used  in  this  work  are  well  known  for  high  accuracy  or  spectral  convergence. 
However,  spectral  convergence  is  achieved  only  when  the  spectral  methods  are  applied  to 
sufficiently  smooth  problems.  The  loss  of  fast  convergence  for  the  problems  that  have 
potential  to  develop  non-smooth  solutions  in  finite  time  undermines  the  advantages  of 
spectral  methods.  Gottlieb  and  Hesthaven  [16]  reported  the  use  of  spectral  filters  acting  in 
a  similar  fashion  as  additional  dissipative  terms  that  help  to  stabilize  simulations  without 
affecting  the  spectral  convergence.  Besides  the  ability  to  retain  the  spectral  convergence 
and  stabilize  the  simulation,  the  filtering  also  decreases  the  aliasing  errors.  According  to 
Gottlieb  and  Hesthaven  [16],  there  is  no  unique  choice  of  the  filter  function  as  long  as 
certain  basic  requirements  arc  satisfied.  The  baseline  DNS  code  employed  for  simulations 
in  this  work  has  been  developed  by  Diamessis  et  al.  [11]  and  employs  the  exponential  filter 


(12) 


where  p  is  the  filter  order,  N  is  the  highest  mode  in  the  spectral  domain,  &/  is  the  filter 
lag  and  a  —  —  loge*/  with  being  the  machine  precision.  The  filter  in  Equation  (12)  is 
intentionally  selected  such  that  it  affects  only  the  high-frequency  modes  and  its  purpose 
is  only  to  stabilize  the  numerical  simulation  and  enhance  the  convergence  rate  of  the 
approximation. 

The  filtered  solution  fF  may  now  be  expressed  in  terms  of  the  modes  in  Legendre  space 
of  the  numerical  solution  as 


Afc-i 


fF  =  £^)M(zy). 


(13) 
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Figure  1.  One-dimensional  transfer  function.  Dash-dot  line:  explicit  primary  filter;  dashed  line: 
secondary  filter  with  N  =  3;  solid  line:  secondary  filter  with  N  =  5;  dotted  line:  spectral  filter  as  in 
Equation  (12). 

where  kj  is  the  jth  discrete  Legendre  mode.  The  filtering  in  Founer  space  is  done  in  the 
similar  fashion  [16]. 

Transfer  functions  for  different  filters  described  above  are  shown  in  Figure  1 .  By  design, 
the  exponential  filter  affects  only  modes  in  the  vicinity  of  the  mesh  cutoff.  The  secondary 
filter  in  Equation  (4)  with  N  =  5  retains  almost  all  information  for  scales  k  <  where 
k *  =  n /{ A)  is  the  nominal  cutoff  wave  number  for  the  physical  LES  scales.  This  is  the  an 
important  characteristic  of  the  secondary  filter  compared  w  ith  other  candidate  filters  and 
is  chosen  on  this  basis  for  the  current  work.  The  scales  between  k^  and  the  mesh  cutoff 
k%  =  n/h  are  strongly  attenuated  by  the  filter  and  play  a  role  of  the  estimated  SGS  in  TNS. 

The  TNS  equations  are  equivalent  to  the  under- resolved  DNS  stabilized  by  the  ex¬ 
ponential  filter  given  in  Equation  (12).  In  the  under-resolved  DNS,  the  energy  begins  to 
accumulate  in  small  scales  and  the  long-time  dynamics  will  be  incorrect.  The  velocity  esti¬ 
mation  model  [3]  removes  this  accumulated  energy  by  periodic  filtering  w  ith  the  secondary 
filter  of  the  form  given  in  Equation  (4).  The  time  interval  between  applications  of  the  filter 
has  to  be  manually  prescribed  by  the  user  in  advance.  The  method  would  gain  in  generality 
if  the  filtering  operation  could  be  activated  automatically  based  on  a  physically  established 
criterion. 

To  establish  such  a  criterion,  the  combinations  of  an  energy  spectrum  E(k)  for  the 
full  turbulent  velocity  field  and  two  secondary  filters  are  theoretically  considered.  The 
secondary  filters  are  built  from  the  primary  filters  with  different  filter  widths.  We  choose 
two  filter  widths,  h  and  A,  where  h  is  the  mesh  size  used  in  the  simulation  and  A  =  2 h. 
The  filter  width  h  is  related  to  the  grid  cutoff  wave  number,  k J,  by  the  relation  k^  =  7T/h, 
and  the  filter  width  A  =  2h  is  related  to  the  LES  cutoff  wave  number,  k by  the  relation 
kf  =  7i /{ A).  The  filtered  quantities  obtained  using  the  filters  corresponding  to  the  cutoff 
wave  numbers  k%  and  k *  are  denoted  by  tilde  and  hat,  respectively.  Figure  2  shows  the 
attenuated  spectra  of  the  inertial,  dissipation  and  ranges.  If  the  velocity  field  is  filtered 
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(a) 
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(b) 

Figure  2.  Energy  spectrum  plot  of  the  infinite  inertial  range.  Thick  solid  line,  energy  spectrum;  thin 
solid  line,  filtered  energy  spectrum  with  spherical  tophat  filter;  dotted  line*  filtered  energy  spectrum 
with  secondary  filter,  N  =  5;  dashed  line:  filtered  energy  spectrum,  N  =  3;  dash-dot:  filtered  energy 
spectrum  with  primary  filter;  vertical  dashed  line:  vertical  dash-dot  line,  k *  (a)  energy  spectrum 

plot  with  the  filtered  energy  spectrum;  filter  with  filter  width  A,  and  (b)  energy  spectrum  plot  with 
the  filtered  energy  spectrum;  filter  with  filter  width  h. 
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Table  1.  Filtered  energy  ratios. 


Energy  spectrum 

Secondary  filter;  1(A) 

N  =  5  1(h) 

Hh) 

/(A) 

Dissipation 

0.0896 

0.0061767 

0.00689363 

Inertial 

0.0234 

0.00021466 

0.009173504 

Batchelor 

0.0699 

0.00068 

0.009728183 

with  the  former  filter,  the  corresponding  spectrum  of  the  filtered  field  is  denoted  as  E{k) 
and  the  energy  removed  by  the  filtering  in  a  range  of  wave  numbers  resolved  in  TNS  up  to 
khc  =  n/h  is 


rkc 

1(h)  =  /  ‘  (E(k)-  E(k))dk, 

Jo 

and  the  similar  integral  for  the  other  filter.  The  expression 

m  (E(k)  -  Eimk 
/o?  (£(*)  -  m)dk 

is  the  ratio  of  energies  removed  by  both  filters  in  the  range  of  wave  numbers  resolved  in 
TNS.  This  ratio  has  been  computed  for  the  progressively  shallower  energy  spectra,  starting 
with  the  dissipation  range  spectrum  form  [17],  the  inertial  range  k~ 5/3  spectrum  and  the 
Batchelor’s  k~[  spectrum.  For  illustration,  in  Figure  2(a)  and  (b)  we  plot  the  inertial  range 
spectrum  and  the  filtered  energy  spectra  for  two  different  filter  widths.  The  results  for  the 
computed  energy  ratio  are  shown  in  Table  1.  Depending  on  the  spectrum,  this  ratio  varies 
in  the  narrow  range  of  between  0,007  and  0.010.  We  will  use  these  bounds  as  guidelines  to 
decide  when  the  filtering  in  the  simulations  should  be  activated 

In  actual  TNS  for  turbulent  channel  flow,  the  energy  ratio  is  computed  from  the  simu¬ 
lation  data  using  the  formula 


(14) 


(15) 


lm\  (")  —  (  2  ^n)]\ 

\I(A)lsiM  \f  £*=:[(«-.  -«nK«n  -  M/.)]  I SIM 

where  {)  is  a  plane  averaged  quantity,  followed  by  integrating  it  over  z ■ 


(16) 


(17) 


When  this  last  quantity  exceeds  the  theoretical  value  (1(h)/  I(A))Theory  from  Table  2,  it  is 
an  indication  that  the  small  scales  in  TNS  become  too  energetic.  Therefore,  when 


m\  ( m  \ 

/(A)/ 

SIM  >  \l( A)J  T hcory 


(18) 


the  filtering  will  be  activated  to  attenuate  the  small  scales.  Thus,  this  condition  serves  as 
the  automatic  filtering  criterion  in  TNS. 
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3.  Numerical  simulations 

The  method  is  validated  in  this  study  by  performing  LES  of  the  incompressible  turbulent 
channel  flow  because  of  the  wealth  of  the  experimental  and  the  DNS  data  that  can  be  used 
for  comparison.  The  fluid  is  contained  between  two  solid  walls  separated  by  the  distance 
Lz.  The  horizontal  domain  has  the  dimension  Lx  for  the  stream  wise  direction  and  L  y  for  the 
spanwise  direction.  The  no-slip  boundary  conditions  are  imposed  on  the  upper  and  lower 
walls  and  penodic  boundary  conditions  are  applied  in  the  horizontal  directions.  The  code 
uses  a  pseudospectral  numerical  method  based  on  the  Fourier  expansions  in  the  streamwise 
and  the  spanwise  directions  and  the  Legendre  polynomials  in  the  vertical  direction  (see 
[11]).  In  the  numerical  method  used  in  this  paper,  the  spectral  filtering  is  applied  in  both 
the  Fourier  and  the  Legendre  spaces  to  maintain  the  stability  and  spectral  accuracy  of  the 
solutions,  and  also  to  help  in  eliminating  the  aliasing  effects.  The  MPI-based  parallel  solver 
was  employed  in  the  numerical  simulations  presented  in  this  paper.  The  assignment  of 
different  sections  of  the  computational  domain  to  individual  processors  is  based  on  one¬ 
dimensional  domain  decomposition,  which  partitions  the  domain  in  distinct  vertical  slabs 
of  thickness  Lx/Np  and  Ly/Np,  when  operating  in  physical  (Fourier)  space,  where  Np  is 
the  number  of  processors.  The  simulations  were  performed  at  the  University  of  Southern 
California  High  Performance  Computing  Center  s  Linux  cluster.  In  order  to  reduce  the  total 
computational  time  and  cost  caused  by  the  interprocessor  communication,  the  condition 
for  the  automatic  filtering  indicated  in  Equation  (18)  was  checked  at  every  20  time  steps, 
while  the  spectral  filtering  with  the  exponential  filter  given  in  Equation  (12)  was  applied  at 
every  time  step  in  the  simulations. 

The  TNS  equations  are  the  standard  Navier-Stokes  equations  for  incompressible  flows 


du{ 

“a7 


duiUj 

dXj 


dUi 

dXf 


1  dp  32Ui 

p  dxt-  dXjdXj  ’ 


0, 


(19) 

(20) 


where  the  tilde  notation  is  used  to  indicate  that  the  spectral  support  for  a  given  quantity  in 
discretized  equations  is  not  sufficient  to  capture  all  dynamically  relevant  scales  of  motion  but 
is  sufficient  to  resolve  fields  filtered  with  the  filter  given  in  Equation  (4).  The  pressure  term  is 
decomposed  into  mean  P(x)  and  fluctuating  component  p(x,  ().  For  a  given,  constant  mean 
pressure  gradient  driving  the  flow,  the  statistically  steady  state  is  eventually  established  in 
which  r0  =  —S(dP/dx)y  where  r0  is  the  wall  shear  stress  and  S  —  Lz/2  is  the  channel  half 
width.  By  introducing  the  shear  velocity  ux  =  ( To/p )^2,  the  Navier-Stokes  equations  are 
rewritten  as  follows: 


dutUj  /  1  dp  u2  $  \  d2w, 
dt  +  dx i  dx{  S  /  dxjdxj 


(21) 


A  statistically  steady  state  is  established  much  faster  if  the  constant  mass  flux  is  imposed 
instead  of  the  constant  pressure  gradient.  The  constant  mass  flux  in  the  simulations  is 
enforced  by  setting  the  mean  pressure  gradient  to  the  instantaneous  wall  shear  stress  at  each 
time  step,  which  is  equivalent  to  using  in  Equation  (21)  the  following  expression  for  uT 


d(«i) 

dx} 


(Lz<  t)+  v 


HSi) 

dX} 


(0,r), 


«r(0  = 


(22) 
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where  (...)  indicates  the  plane  averaged  quantity.  The  above  equations  are  normally  nondi- 
mensionalized  by  the  channel  half  width  8  =  Lzj 2  and  the  nominal  friction  velocity 
«o  =  (ro//?)l/2  for  the  constant  pressure  gradient  case,  i.e.  Tq  =  —  8(dP/dx),  giving 


duj  dujUj  _  dp_  \  ,  1 

St  dxi  'p3x/  T  lJ )  Re0  3xj3xj 1 


(23) 

(24) 


where  Re o  =  ^  is  the  nominal  Reynolds  number.  Note  that  the  nominal  Reynolds  number 
Re o  is  used  only  as  the  parameter  in  the  simulations  and  its  value,  while  close,  is  usually 
different  from  the  actual  Reynolds  number  Rex  =  ux8/v  based  on  the  actual  friction  veloc¬ 
ity  as  shown  in  Equation  (22).  Only  in  the  statistically  steady  state  and  in  constant  pressure 
gradient  simulations  Re o  =  Rex .  The  actual  friction  velocity  nondimensionalized  by  uq  is 


«r(0 


l  Hu i) 
Re o  dx 3 


(2.0  + 


1  d(u  |) 


(0,0- 


(25) 


The  parameters  used  in  the  simulations  are  gathered  in  Table  2  and  the  cases  used  for 
comparisons  are  collected  in  Table  3.  The  LES  simulations  were  run  until  they  reached  the 
statistically  steady  state.  In  each  case,  the  filtering  criterion  has  been  selected  from  the  the 
range  between  0.007  and  0.009  suggested  by  the  theory.  The  comparison  data  in  Table  3 
are  taken  from  the  literature.  Cases  TNSZ1  and  TNSZ2  are  the  results  from  the  LES  using 
the  velocity  estimation  model  by  Domaradzki  et  al.  [3].  The  case  HiDNS  is  from  Gilbert 
[18]  at  Rex  =210.  The  DNS  results  at  Rex  =  2003  for  the  case  HiDNS2  are  from  the  work 
of  Hoyas  et  al.  [19].  The  case  HiDNS3  are  the  DNS  results  at  Rex  —  944  from  the  work 
of  Del  Alamo  et  al.  [20].  The  ratio  of  the  total  number  of  mesh  points  used  in  TNS  and 
DNS  shown  in  Table  2  demonstrates  that  TNS  resolutions  can  be  of  the  order  of  few  tenths 
of  1%  of  the  DNS  resolutions.  This  indieates  significant  saving  in  term  of  computational 
effort  for  TNS  compared  with  DNS. 


4.  Results  and  discussion 

In  this  section,  the  results  from  the  simulations  for  all  cases  from  Table  2  are  compared 
with  the  baseline  data  for  cases  listed  in  Table  3.  The  comparisons  involve  mean  and  RMS 
turbulent  velocities. 

The  Re  180  case  is  the  classical  test  case  at  Rex  ~  180  w  ith  the  domain  size  4n  x  2n  x  2 
for  which  the  detailed  DNS  data  exist  [18].  The  filtering  condition  is  set  to  0.008  resulting 
in  the  filter  being  turned  on  automatically  every  360  to  580  time  steps.  The  case  is  run 
for  24,000  times  steps  and  is  well  converged  to  a  statistically  steady  state  when  the  data 
for  plotting  are  taken.  In  Figure  3(a)  we  plot  the  mean  velocity,  and  the  RMS  velocities 
are  shown  in  Figures  3(b)— (d).  The  results  for  the  case  Re  180  are  compared  with  the 
high  resolution  DNS  results  of  Gilbert  and  Kleiser  [18]  (case  HiDNS),  TNS  results  of 
Domaradzki  et  al.  [3]  (case  TNSZ1)  and  under-resolved  DNS  performed  with  the  same 
resolution  of  64  x  64  x  65  as  for  the  LES  cases  Re  180  and  TNSZ1.  The  mean  velocity  in 
the  Re  180  case  agrees  quite  well  with  both  HiDNS  and  TNSZl  cases.  Only  in  the  range  of 
Z +  between  7  and  1 5,  the  result  from  this  simulation  slightly  over  predicts  results  from  those 
two  baseline  cases.  This  overshoot  of  the  mean  velocity  is  due  to  the  slight  under-prediction 
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Table  3  Parameters  for  the  reference  simulations. 


Case 

Domain  size 

Resolution 

Ajc+  x  A y+  x  A z+ 

Re, 

TNSZ1 

4tt  x  2n  x  2 

64  x  64  x  65 

35.73  x  17.88  x  0  85 

182 

TNSZ2 

2.5 7T  X  0.57T  x  2 

96  x  128  x  65 

78.62  x  11.79  x  4.49 

961 

HiDNS 

3.6;r  x  1.97T  x  2 

160  x  160  x  129 

14.85  x  7.83  x  0.063 

210 

HiDNS2 

87r  x  37T  x  2 

6144  x  4608  x  633 

8.3  x  44  x  8.9a 

2003 

HiDNS3 

87t  x  3tt  x  2 

3072  x  2304  x  385 

7.6  x  7.6  x  0.032 

944 

Note:  aFor  HiDNS2,  Az+  =  Az*a;r,  for  all  olher  cases  Az+  is  lhe  distance  between  the  first  grid 
point  and  the  boundary 


of  the  energy  dissipation  in  the  buffer  region.  The  imbalance  between  the  kinetic  energy 
production  and  the  dissipation  is  caused  by  less  frequent  filtering  compared  with  the  TNSZ 1 
case.  The  under- resolved  DNS  underpredicts  the  mean  velocity  in  the  log  region.  In  Figure 
3(b),  the  streamwise  RMS  velocity  in  the  simulation  is  plotted  together  with  all  companson 
cases  as  in  Figure  3(a).  The  urms  from  simulation  compares  well  with  both  baseline  cases 
HiDNS  and  TNSZ1.  The  same  trend  continues  for  the  spanwise  RMS  velocity  shown  in 
Figure  3(c),  but  the  vertical  RMS  velocity  is  somewhat  underpredicted  for  both  Re  180  and 
TNSZ1  cases  (Figure  3(d)). 

The  second  test  case  Re  1050  is  LES  performed  in  a  domain  w  ith  the  size  2.5;r  x  0.5;r  x 
2,  the  resolution  of  96  x  128  x  65  mesh  points,  at  Re o  =  1050.  The  results  for  this  case  are 
plotted  in  Figure  4  together  with  the  experimental  results  from  Wei  and  Wilmarth  [21],  the 
DNS  results  from  Del  Alamo  et  al.  [20],  the  LES  results  from  Piomelli  [22]  obtained  using 
the  dynamic  Smagorinsky  model,  the  LES  results  from  Domaradzki  et  al.  [3]  obtained 
with  the  velocity  estimation  model  (the  case  TNSZ2)  and  the  results  from  DNS  without 
any  modeling  but  at  the  same  low  resolution  of  96  x  128  x  65  as  in  two  LES  cases.  The 
resolution  chosen  is  the  same  as  the  one  used  for  the  case  TNSZ2  in  Domaradzki  et  al.  [3] 
and  comparable  to  the  resolution  used  by  Piomel  li  [22].  The  filtering  condition  is  set  to  0.007 
and  the  filter  is  found  to  be  activated  every  350  to  620  time  steps.  The  statistically  steady 
state  was  reached  after  about  16,000  time  steps.  All  LES  mean  velocity  data,  including 
the  case  Re  1050,  are  shown  in  Figure  4(a)  and  are  in  good  mutual  agreement,  though  they 
slightly  overpredict  the  experimental  data.  Figure  4(b)  shows  the  plot  of  the  streamwise 
RMS  velocity  from  simulation  compared  with  other  baseline  cases.  The  agreement  with  the 
experimental  data  and  the  TNSZ2  data  is  good  between  z  =  0  andz  =  0.2,  but  the  quality  of 
comparison  deteriorates  for  larger  values  ofz.  The  vertical  RMS  velocity  LES  data  shown 
in  Figure  4(c)  are  in  good  agreement  with  the  experiments  throughout  the  entire  domain. 
No  results  for  the  spanwise  velocity  are  presented  here  because  the  experimental  data  for 
that  velocity  were  not  available  for  comparisons. 

This  Reynolds  number  was  also  selected  to  investigate  in  more  detail  the  dependence 
of  TNS  on  the  mesh  resolution.  It  is  well  known  (e.g.  [22])  that  the  predictions  of  the 
mean  velocity  profile,  RMS  velocities  and  high-order  statistics  depend  on  the  gnd  resolu¬ 
tion,  particularly  in  the  wall  region.  We  have  performed  two  additional  TNS  simulations 
at  Re o  =  1050  with  increased  resolutions.  The  case  Rel050h  is  the  LES  performed  in  the 
same  domain  size  as  the  one  used  in  the  case  Rel050,  but  with  the  grid  resolution  adjusted 
according  to  the  discussion  in  Piomelli  [23],  Piomelli  and  Balaras  [24]  and  the  recommen¬ 
dations  by  Sagaut  [25].  The  main  effect  of  the  changed  resolution  is  that  the  first  mesh  point 
away  from  the  wall  is  at  z+  <  1 ,  as  required  for  the  wall-resolved  LES.  The  Figures  5(a)-(c) 
show  the  plots  obtained  from  the  simulation  for  this  case  compared  with  the  experimental 
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and  DNS  data  from  the  same  sources  as  the  reference  data  used  in  Figures  4(a)-(c).  Thus, 
the  mean  velocity  profile  shown  in  Figure  5(a)  improves.  Also,  for  this  resolution,  the  mean 
velocity  for  the  under- resolved  case  is  predicted  quite  well.  This  case  uses  only  the  weak 
stability  filter  of  Equation  (12)  at  each  time  step  and,  based  on  the  prediction  of  the  mean 
velocity,  provides  a  sufficient  substitute  for  an  explicit  SGS  model.  The  TNS  filter  leads  to 
noticeable  improvements  in  the  prediction  of  the  RMS  velocities  over  the  lower  resolution 


Figure  3.  Results  from  simulation  of  the  case  Re  180.  Dashed  line:  Re  180;  diamond  mark*  HiDNS; 
dotted  line:  TNSZ1;  dash-dot  line:  under-resolved  DNS;  (a)  mean  velocity  profiles,  (b)  streamwise 
RMS  turbulent  velocity,  (c)  spanwisc  RMS  turbulent  velocity  and  (d)  vertical  RMS  turbulent  velocity. 
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Figure  3.  (Continued) 

case  Rcl050.  The  results  obtained  using  the  TNS  filtering  arc  also  better  than  For  the  under- 
resolved  case,  though  not  much.  The  effective  friction  Reynolds  number  also  shows  the 
improvement  in  capturing  the  near- wall  dynamics  by  the  computational  grid.  However,  there 
is  still  the  overshoot  of  the  streamwise  component  of  the  RMS  velocity  in  near-wall  region 
and  the  wall  normal  component  of  the  RMS  velocity  in  near-wall  region  is  slightly  under¬ 
predicted,  indicating  that  the  gird  does  not  fully  resolve  the  structures  in  the  near-wall  region 
which  are  significant  in  the  turbulent  energy  production  and  for  the  SGS  interactions.  The 
case  Rel050h2  doubles  the  mesh  resolution  in  the  streamwise  direction  to  satisfy  the  gener¬ 
ally  accepted  criteria  for  wall-rcsolvcd  LES  [23-25].  The  results  arc  significantly  improved 
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for  the  TNS  case,  but  no  further  improvement  is  observed  for  the  under- resolved  case  with 
the  stability  filter.  The  mean  velocity  profile  and  the  RMS  velocity  plots  agree  very  well  with 
the  DNS  data  as  seen  in  the  Figures  6(a)-(c).  The  small  overshoot  in  the  buffer  region  is  still 
present  but  the  turbulent  intensities  in  the  log  region  are  different  from  the  intensities  in  the 
cases  Re  1050  and  Rel050h,  There  is  no  overshoot  in  near- wall  region  for  RMS  streamwise 
component  of  the  velocity  and  the  RMS  wall-normal  component  of  the  velocity  agrees  very 
well  with  the  DNS  data.  These  results  indicate  that  the  grid  resolution  in  the  Rel050h2 
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Figure  4.  Results  from  simulation  of  the  case  Re  1050.  Dashed  line:  Re  1050,  triangle  mark*  Wei  & 
Wilmart;  dash-dot  line  with  plus  mark:  Piomelh;  dotted  line:  TNSZ2,  dash-dot  line:  under-resolved 
DNS;  diamond  mark:  HiDNS3;  (a)  mean  velocity  profiles,  (b)  streamwise  RMS  turbulent  velocity 
and  (c)  vertical  RMS  turbulent  velocity 
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Figure  4.  (Continued) 


case  in  conjunction  with  the  TNS  procedure  is  able  to  resolve  the  near-wall  dynamics. 
However,  though  the  application  of  only  stability  filter  allows  to  capture  the  mean  flow,  the 
prediction  of  the  rms  velocities  for  this  case  does  not  have  the  same  quality  as  the  full  TNS 
method. 

Figures  7(a)~(d)  show  plots  obtained  in  LES  at  Re o  =  2000  in  the  domain  size  2.5 n  x 
0.5tt  x  2,  and  the  resolution  256  x  128  x  131  The  filtering  criterion  is  set  to  0  008  and 
the  filter  turns  on  automatically  every  80  to  160  time  steps.  This  is  more  frequent  filtering 
than  in  previous  cases,  indicating  that  the  simulations  at  high  Reynolds  number  require 
more  dissipation.  The  statistically  steady  state  is  reached  after  about  21, 000  time  steps. 
The  comparisons  are  made  with  the  under-resolved  DNS  and  the  HiDNS2  case,  which  is 
the  fully  resolved  DNS  of  Hoyas  and  Jimenez  [19],  but  in  the  larger  domain  87T  x  3;r  x  2. 
The  resolution  used  in  the  Re2000Sh  case  is  selected  based  on  the  observations  from  the 
cases  Rel050h  and  Rel050h2  and  also  from  the  computational  expense  considerations. 
The  mean  velocity  shown  in  Figure  7(a)  slightly  overpredicts  the  DNS  results  in  the 
similar  way  as  in  the  lower  Reynolds  number  case  Rel050h  (Figure  5(a)).  The  quality  of 
the  prediction  for  the  streamwise  RMS  velocity  presented  in  Figure  7(b)  is  lower,  LES 
overpredicts  the  peak  of  the  urms  from  the  HiDNS2  case  and  underpredicts  the  DNS  results 
in  the  range  >  200.  The  remaining  two  RMS  turbulent  velocity  components  are  shown  in 
Figures  7(c)  and  exhibit  good  agreement  with  the  DNS  data  at  z+  <  800.  The  overprediction 
and  the  underprediction  of  the  peak  of  urms  could  come  from  the  insufficient  grid  resolutions 
in  the  near-wall  region.  Nevertheless,  at  this  high  Reynolds  number,  for  the  computational 
expense  the  quality  of  the  low-order  statistics  is  quite  satisfactory. 
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Figure  5.  Results  from  simulation  of  the  case  Rel050h.  Dashed  line:  Rel  050h;  triangle  mark:  Wei 
&  Wilmart;  dash-dot  line  with  plus  mark:  Piomclli;  dotted  line:  TNSZ2;  dash-dot  line:  under- resolved 
DNS;  diamond  mark.  HiDNS3;  (a)  mean  velocity  profiles,  (b)  streamwise  RMS  turbulent  velocity 
and  (c)  vertical  RMS  turbulent  velocity. 
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Figure  6.  Results  from  simulation  of  the  case  Re  1 050h2.  Dashed  line:  Re  1 050h2;  triangle  mark:  Wei 
&  Wilmart;  dash-dot  line  with  plus  mark:  Piomelli;  dotted  line:  TNSZ2;  dash-dot  line:  under-resolved 
DNS;  diamond  mark:  HiDNS3;  (a)  mean  velocity  profiles,  (b)  streamwise  RMS  turbulent  velocity 
and  (c)  vertical  RMS  turbulent  velocity. 
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Figure  6.  (Continued) 
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(b) 

Figure  7,  Results  from  simulation  of  the  case  Re2000Sh,  Dashed  line:  Re2000Sh;  diamond  mark: 
HiDNS2;  dash-dot  line:  Under  resolved  DNS  (a)  Mean  velocity  profiles,  (b)  Streamwise  rms  turbulent 
velocity,  (c)  Spanwise  rms  turbulent  velocity,  and  (d)  Vertical  rms  turbulent  velocity 
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5.  Conclusions 

The  numerical  code  employed  in  this  work  uses  spectral  methods  stabilized  by  a  weak  filter 
which  allows  it  to  run  for  all  Reynolds  numbers  and  resolutions  considered.  Despite  the  sim¬ 
ulations  being  stable,  the  results  from  such  runs  do  not  compare  well  with  the  baseline  DNS 
and  the  experimental  data,  as  expected  from  under-resolved  DNS,  unless  very  high  resolu¬ 
tion  is  employed.  It  is  thus  clear  that,  beyond  numerical  stabilization  techniques,  additional 
turbulence  modeling  procedures  are  required  to  obtain  physically  correct  results  in  simu¬ 
lations  that  do  not  resolve  all  relevant  scales  of  motion.  We  have  presented  the  LES  model 
that  utilizes  the  low-pass  filtering  operation  as  the  only  modeling  tool.  The  entire  procedure 
follows  closely  the  previously  developed  TNS  approach  which  uses  the  filter  periodically 
to  modify  smallest  resolved  scales.  In  the  present  TNS  implementation,  the  energy  of  these 
scales  is  monitored  and  the  filter  is  turned  on  automatically  whenever  the  energy  exceeds  a 
threshold  defined  by  a  specific  criterion.  We  have  tested  the  model  on  the  turbulent  channel 
flow  problem  at  different  Reynolds  numbers,  computational  resolutions  and  domain  sizes. 
Satisfactory  comparisons  among  TNS  results  and  various  reference  cases  have  been  ob¬ 
tained.  In  particular,  the  quality  of  the  results  is  comparable  to  the  previous  TNS  procedure 
which  uses  a  fixed  time  interval  for  the  filtering  as  well  as  additional  SGS-estimated  velocity. 

The  main  goal  of  LES  is  to  predict  accurately  low-order  statistics  using  the  lowest  possi¬ 
ble  numerical  resolution.  Because  of  that  the  grid  resolutions  used  in  the  reported  work  were 
selected  to  test  and  study  the  dependency  of  the  quality  of  low-order  statistics'  predictions  on 
the  grid  resolutions,  using  guidelines  from  previous  work  regarding  the  wall-resolved  LES 
of  the  channel  flow.  The  cases  Re  1 050  and  Re2000Sh  were  selected  to  test  the  proposed  TN  S 
method  with  coarse  resolution.  The  overpredicted  wall  shear  stress  and  the  overshoot  in  the 
near-wall  region  in  urms  show  that  the  grid  resolutions  are  not  sufficient  to  resolve  the  struc¬ 
tures  near  the  wall.  The  TNS  filtering  that  provides  SGS  dissipation  tends  to  overdissipate  en¬ 
ergy  in  the  log  region  as  seen  in  the  results  forwrmi  and  wrms  for  these  cases.  When  resolution 
is  increased,  the  quality  of  the  results  improves  significantly.  This  is  clearly  seen  in  the  plots 
of  Urms  for  foe  cases  Re  1050,  Rel050h  and  R1050h2,  where  the  overshoot  in  the  near-wall 
region  decreases  when  the  grid  resolution  increases.  The  mean  velocity  profiles  for  these  test 
cases  also  improve  and  this  improvement  makes  the  shear  stress  predictions  more  accurate. 
The  case  Re2000Sh  was  conducted  to  test  the  performance  of  our  approach  at  high  Reynolds 
number  as  the  simulations  would  expect  to  depend  more  on  the  filtering  to  provide  energy 
dissipation.  Even  though  the  simulations  in  these  cases  are  done  at  high  Reynolds  number 
and  the  employed  gnd  resolutions  arc  marginal,  the  predictions  of  the  mean  velocity  profiles 
are  reasonable.  The  improvement  in  the  quality  of  the  predictions  using  this  approach  at 
high  Reynolds  numbers  should  be  expected  as  the  grid  resolution  is  increased  further. 

In  summary,  we  have  shown  that  the  TNS  method  is  capable  to  produce  quality  LES 
results  despite  a  lack  of  an  explicit  SGS  model.  In  TNS,  the  smallest  scales  resolved  by 
the  mesh  provide  a  model  of  SGS.  They  are  modified  periodically  by  the  explicit  filtering 
operation,  which  is  triggered  by  a  condition  based  on  the  analysis  of  the  energy  spectra 
in  developed  turbulence.  We  believe  that  the  simplicity  and  the  flexibility  are  the  features 
of  the  method  that  should  be  appreciated.  Like  with  any  modeling  procedure,  there  is 
always  room  for  improvement.  For  instance,  it  is  rather  obvious  that  the  filtering  criterion 
will  depend  on  the  filter  s  form.  However,  one  may  conjecture  that  TNS  results  should  be 
independent  of  the  filter’s  form  as  long  as  the  filter’s  effect  is  to  substantially  attenuate  small 
scales  while  leaving  the  large  scales  unchanged  and  if  the  filtering  criterion  is  derived  for 
the  specific  filter  used.  We  will  investigate  the  dependence  of  the  model  on  the  filter  form 
in  the  subsequent  work. 
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1.  Introduction 

The  equations  for  large  eddy  simulation  (LES)  are  traditionally  obtained  by  the  filtering 
approach  proposed  by  Leonard  (1974)  where  a  spatial  filter  that  strongly  attenuates 
scales  of  motion  smaller  than  the  prescribed  filter  width  A  is  applied  to  the  Navier- 
Stokcs  equations.  If  the  filtered  quantities  arc  denoted  by  an  overbar,  the  LES  equations 
for  an  incompressible  flow  become 


0  _  d _  Id  d2  __  d 

■^ut  +  -X—  Ui  Uj  =  — —p  +  v- — —Ui  ~  —  rtJ, 
at  axj  paxx  dXjOXj  ox3 


a) 


dnUi  ~  0  ’  (2) 

where  ux  =  (ui,  U2,u$)  =  (u,vxw)x  p,  and  v  are  the  velocity,  pressure,  and  the  kinematic 
viscosity,  respectively,  and  rl3  is  the  subgrid-scale  (SGS)  stress  tensor 

Tij  =  TT[u3  —  Ui  u3 .  (3) 

The  form  of  equations  (1)  and  (2)  requires  that  the  filtering  and  differentiation  com¬ 
mute  (Ghosal  and  Moin  (1995),  Vasilyev  et  ai  (1998)).  In  practice,  however,  the  above 
equations  are  frequently  the  starting  point  in  SGS  modeling  without  regard  to  formal 
requirements  for  their  derivation  in  the  filtering  framework.  The  important  point  is  that 
the  LES  equations  have  the  form  of  the  Navier-Stokes  equations  for  the  filtered  velocity 
ux  plus  the  additional  force  term  which  is  the  divergence  of  the  subgrid  scale  stress  tensor 
(3),  and  which  is  required  to  close  the  LES  equations.  Various  SGS  models  differ  in  how 
the  SGS  stress  tensor  is  expressed  (or  modeled)  in  terms  of  the  filtered  velocity  ut . 

Among  SGS  models  the  most  important  category  are  the  eddy  viscosity  models.  Their 
origin  goes  back  to  Boussinesq  (1877)  who  proposed  that  the  effects  of  turbulence  can 
be  accounted  for  by  the  viscosity  increased  over  its  value  in  laminar  flows.  With  the 
advent  of  computers  and  attempts  to  perform  general  circulation  simulations,  the  eddy 
viscosity  concept  was  used  by  Smagorinsky  (1963)  to  model  uThe  lateral  transfer  of 
momentum  and  heat  by  the  non-linear  diffusion,  which  parametrically  is  supposed  to 
simulate  the  action  of  motions  of  sub-grid  scale,...”,  thus  starting  the  modern  era  of  LES 
and  SGS  modeling.  There  are  a  number  of  excellent  reviews  of  theory  and  practice  of 
SGS  modeling  in  the  eddy  viscosity  framework,  e.g.,  by  Galperin  and  Orszag  (1993), 
Lesieur  and  Metais  (1996),  Piomelii  (1999),  Meneveau  and  Katz  (2000),  Pope  (2000), 
and  Sagaut  (2002),  and  such  methods  are  not  addressed  in  this  review. 

Despite  its  practical  utility  the  eddy  viscosity  concept  has  not  been  accepted  as  a  fully 
satisfactory  solution  of  the  turbulence  modeling  problem  because  it  is  not  consistent  with 
some  observed  features  of  turbulence  physics.  For  instance,  a  purely  dissipative  character 
of  the  eddy  viscosity  does  not  allow  such  models  to  describe  an  inverse  energy  trans¬ 
fer  (backscatter)  observed  in  actual  turbulent  flows.  The  dynamic  Smagorinsky  model 
of  Germano  et  al.  (1991)  formally  allows  for  the  backscatter  but  this  may  lead  to  unstable 
simulations.  In  practice  various  averaging  procedures  for  the  dynamic  Smagorinsky  coef¬ 
ficient  are  used  to  render  the  eddy  viscosity  a  positive  quantity.  Nevertheless  the  dynamic 
procedure  offers  a  significant  improvement  over  the  eddy  viscosity  models  with  constant 
coefficients.  This  is  because  the  modeled  SGS  transfer  is  consistent  with  the  actual  en- 
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ergy  transfer  across  the  test  cutoff,  allowing  the  model  to  better  represent  energy  drain 
from  the  resolved  scales  for  a  wider  variety  of  turbulent  flows.  Yet  this  desirable  feature 
may  turn  into  disadvantage  in  some  situations.  One  such  situation  is  turbulence  with 
rotation.  Traditional  eddy  viscosity  models  encounter  a  number  of  difficulties  when  ap¬ 
plied  to  rotating  flows  (Speziale  (1989),  Squires  and  Piomelli  (1995))  because  the  energy 
transfer  from  the  large  scales  to  the  small  ones  is  reduced,  and  consequently  the  energy 
dissipation  decreases  with  rotation  rate  (e.g.,  Morinishi  et  ai  (2001),  Yeung  and  Zhou 
(1998),  Jacquin  et  ai  (1990)),  requiring  lower  values  of  the  Smagorinsky  constant  than 
for  non-rotating  flows.  The  dynamic  procedure  accounts  for  this  effect  but  it  also  makes 
the  Smagorinsky  coefficient  a  spatially  variable  quantity.  As  shown  by  Horiuti  (2001) 
and  Kobayashi  and  Shimomura  (2001)  because  of  the  latter  feature  the  model  violates 
transformation  properties  between  inertial  and  rotating  frame  of  reference.  Rather  than 
using  such  specific  arguments  sometimes  more  general  arguments  are  invoked  against 
eddy  viscosity  models.  For  instance,  one  may  note  that  LES  attempts  to  simulate  scales 
which,  by  virtue  of  being  far  removed  from  the  dissipation  range,  are  essentially  inviscid; 
simulating  them  as  dissipative  scales  through  the  use  of  an  eddy  viscosity  appears  as  the 
contradiction.  Listing  such  few  limitations  and  contradictions  is  not  meant  to  disqualify 
the  eddy  viscosity  models.  Indeed,  the  experience  accumulated  over  many  decades  in 
this  field  is  that  the  single  most  important  requirement  in  turbulence  modeling  is  that  a 
model  is  globally  dissipative,  and  the  eddy  viscosity  concept  guarantees  this  by  design. 
However,  the  recognition  of  the  above  listed  difficulties  serves  to  motivate  research  on 
other  approaches  to  the  SGS  modeling  that  seek  improvements  in  the  fidelity  of  LES  pre¬ 
dictions  without  constraints  imposed  by  the  eddy  viscosity  concept.  The  main  subject 
of  this  paper  is  a  review  of  a  collection  of  approaches  to  the  problem  of  SGS  model¬ 
ing,  which  do  not  use  explicit  eddy  viscosity  expressions.  While  no  single  approach  has 
emerged  to  displace  the  current  eddy  viscosity  models,  the  continuing  progress  in  this 
area  offers  hope  that  such  approaches  may  eventually  become  successful  in  simplifying 
the  task  of  LES  and  in  improving  its  predictive  capabilities. 

2.  Similarity  Models 

In  well  resolved  direct  numerical  simulations  full  velocity  fields  are  available.  For  a  pre¬ 
scribed  filter  the  filtered  velocity  can  be  computed  numerically  from  the  filter  definition, 
often  represented  as  a  convolution  integral  with  a  given  smoothing  kernel  G, 


(4) 


Similarly,  for  a  given  u,  the  exact  SGS  stress  tensor  can  be  computed  from  its  defini¬ 
tion  (3).  Independently,  the  modeled  SGS  tensor  can  be  found  using  the  known  U{  and 
parameters  of  the  model  and  both  stresses  can  be  compared,  constituting  the  so-called 
a  priori  comparison.  It  was  realized  early  by  Clark  et  al.  (1979),  and  later  confirmed  in 
independent  investigations  (Piomelli  et  al.  1988,  Doinaradzki  et  al.  1993,  Liu  et  ai  1994), 
that  the  eddy  viscosity  SGS  stress  expressions  correlate  poorly  with  the  actual  stresses. 
This  conclusion  provided  an  original  impetus  for  the  search  for  the  models  that  would 
improve  the  quality  of  a  prion  comparisons.  The  first  model  that  offered  a  considerable 
improvement  in  this  area  has  been  proposed  by  Bardina  et  al.  (1983)  and  is  known  as 
the  similarity  model.  While  the  original  reasoning  leading  to  the  model  was  more  com¬ 
plicated,  the  final  result  can  be  viewed  as  the  result  of  approximating  the  full,  unknown 
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in  LES  velocity  Ui,  by  the  filtered  velocity  ux  which  is  known  in  LES 


(5) 


and  using  it  in  the  definition  (3),  i.e. 


Tij  =  UxUj  —  Ui  Uj 


ux  Uj  . 


(6) 


The  similarity  model  proposed  by  Liu  et  al.  (1994)  and  O'Neil  and  Meneveau  (1997) 
also  uses  the  known,  resolved  velocity  ux  as  an  approximation  to  the  full  velocity  but 
the  modeled  SGS  stress  tensor  (3)  is  assumed  to  be  proportional  to  the  resolved  stress 
tensor  LXJ  computed  using  a  different  filter,  denoted  here  by  a  hat 


(7) 


The  hat  filter  is  wider  than  the  overbar  filter,  therefore  the  expression  for  Lij  reflects 
qualitatively  the  property  of  the  exact  SGS  stress  (3)  where  the  underlying  velocity  ux 
contains  smaller  scales  than  those  that  are  retained  by  the  application  of  the  filter.  Both 
expressions,  (6)  and  (7),  show  much  improved  correlations  in  numerical  and  experimental 
a  prion  tests. 

However,  the  practical  experience  is  that  when  such  models  are  implemented  in  actual 
LES  the  simulation  results  quickly  deteriorate,  with  the  main  culprit  identified  as  an 
insufficient  SGS  dissipation  (Liu  et  al.  1994).  Because  of  that  a  number  of  so-called 
mixed  models  have  been  proposed  that  simply  combine  similarity  models  (to  guarantee 
a  good  performance  in  a  priori  tests)  and  eddy  viscosity  models  (to  provide  a  good 
performance  in  a  posteriori  tests,  i.e.,  in  actual  LES).  The  solution  proposed  by  the 
mixed  models  is  not  entirely  satisfactory.  It  is  obvious  that  it  is  the  presence  of  the  eddy 
viscosity  expressions  that  makes  the  mixed  models  wrork  in  actual  LES.  Since  comparable 
LES  results  can  be  obtained  with  only  eddy  viscosity  models  it  is  not  clear  that  the 
similarity  component  plays  any  useful  role  in  LES  performed  with  the  mixed  models. 
The  stronger  statement  is  that  the  need  for  adding  the  eddy  viscosity  expression  to  a 
similarity  expression  signifies  the  failure  of  the  similarity  concept  as  a  SGS  modeling 
approach.  The  current  LES  practice  implicitly  endorses  this  view  through  a  diminishing 
emphasis  on  similarity  models  and  a  prion  tests.  Nevertheless,  the  understanding  of  the 
failure  of  the  similarity  modeling  is  useful  in  developing  non-eddy  viscosity  models,  to 
avoid  similar  mistakes  on  the  one  hand,  and  on  the  other  to  decide  if  the  similarity 
concept  is  fundamentally  flawed  or  if  only  its  specific  implementations  are  at  fault. 

The  reasons  for  the  failure  of  the  similarity  models  can  be  qualitatively  understood 
if  the  general  formula  (7)  is  considered  for  sharp  spectral  filters  with  two  wavenumber 
cutoffs  k  >  k ,  where  k  is  also  the  numerical  inesh  cutoff.  This  way  all  active  LES  modes 
are  split  into  two  bands:  band  1  for  0  <  k  <  k  and  band  2  for  k  <  k  <  k.  The  physical 
effect  of  the  expression  u3  -  ux  u3  j  has  been  elucidated  in  investigations  dealing  with 

the  energy  transfer  in  turbulent  flows  (Domaradzki  et  al.  1993,  Kerr  and  Doinaradzki 
1996).  The  first  term  in  the  difference  describes  the  effect  of  nonlinear  interactions  among 
all  resolved  LES  modes  (band  1  and  2)  on  modes  in  the  band  1  (the  widehat  filter).  The 
second  term  describes  the  effect  of  nonlinear  interactions  among  modes  from  the  band  1 
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on  modes  in  band  1  and  2.  Sometimes  it  is  convenient  to  use  the  notation  T7*7™  to  indicate 
the  energy  transfer  to  band  m  caused  by  the  nonlinear  interactions  between  bands  p  and 
q  (Domaradzki  and  Rogallo  1990,  Domaradzki  et  al.  1994).  With  this  notation  the  effect 
of  the  similarity  expression  on  the  energetics  of  the  flow  is  represented  symbolically  as 


(T in  +  T121  +  T221)  -  (Tlu  +  T112)  .  (8) 

In  the  above  expression  the  effects  of  energy  redistribution  within  band  1  (term  Txxx) 
cancel.  The  remaining  terms  describe  the  energy  transfer  from  the  band  1  due  to  inter¬ 
actions  between  band  1  and  2  (T121)  and  interactions  of  modes  in  band  2  (T*221);  they 
result  in  the  energy  loss  by  the  band  1.  The  last  term  T112  is  normally  the  energy  gain 
in  band  2  due  to  interactions  among  modes  in  band  1;  but  because  of  the  negative  sign 
it  is  an  effective  loss  term  for  the  band  2.  The  interscale  energy  transfer  caused  by  the 
complete  nonlinear  term  in  LES  equations,  —d/dx3(TIi  Tij)y  accounts  for  all  combinations 
of  interacting  bands  and  can  be  represented  as  follows 


_|_  rp>12l  _|_  /yi221^  _|_  ^^112  _|_  rj,  122  rp 222^ 

Based  on  observed  magnitudes  and  signs  of  all  terms  in  that  expression  (Domaradzki 
et  al.  1994)  it  has  an  interpretation  consistent  with  the  classical  picture  of  the  turbulent 
energy  cascade:  the  energy  is  removed  from  band  1  (terms  TX2X  and  T221)  and  deposited 
into  band  2  (terms  T112  and  T122)  while  T111  and  T 222  have  redistributive  character. 
The  principal  role  of  a  SGS  model  is  to  counteract  the  accumulation  of  the  energy  in  the 
band  2  caused  by  this  energy  transfer  process,  i.e.,  to  counteract  the  effect  of  terms  T112 
and  Tx 22 .  However,  the  similarity  expression  (8)  counteracts  only  one  of  the  two  active 
terms,  T112,  and  also  removes  the  energy  from  band  1  through  the  action  of  T121  and 
T221,  essentially  doubling  the  energy  transfer  from  band  1  caused  by  the  nonlinear  term 
in  LES  equations.  The  result  is  that  the  model  is  deficient  in  removing  the  energy  from 
the  vicinity  of  the  LES  cutoff  and  unnecessarily  increases  the  energy  loss  by  the  large 
scales  in  band  1  beyond  what  is  already  provided  by  the  resolved  nonlinear  interactions. 
The  quantitative  analysis  can  be  more  complicated  for  different  similarity  expressions 
and  because  of  the  presence  of  the  viscous  and  non-st at ionarity  effects.  Nevertheless, 
the  above  energy  transfer  analysis  provides  the  qualitative  explanation  of  the  failure 
of  the  similarity  concept  in  SGS  modeling.  Such  a  qualitative  understanding  may  offer 
suggestions  for  improvements  in  this  category  of  SGS  models. 


3.  Deconvolution  and  Nonlinear  Models 

The  deconvolution  models  are  based  on  the  observation  that,  under  certain  conditions, 
the  filtering  operation  can  be  inverted,  thus  providing  the  unfiltered  quantity  Ui  from 
the  filtered  one  Tq.  The  convolution  operation  (4)  in  the  Fourier  space  is 


m(k)  =  G(k)u{(k ), 


(10) 


where  for  simplicity  we  assume  a  1-D  filter  and  the  Fourier  transform  is  indicated  by  the 
dependence  on  the  wavenumber  k.  The  inverse  operation  can  be  formally  defined  as 
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iii(k)  =  Ui(k)/G{k),  (11) 

and  is  possible  if  the  division  can  be  performed  and  the  result  is  an  integrable  function. 

In  the  context  of  SGS  modeling  the  inversion  of  filtering  lias  been  proposed  first  by 
Leonard  (1974)  and  later  by  Germano  (1986),  though  not  in  the  form  described  above. 
Leonard  (1974)  showed  that  the  SGS  stress  tensor  can  be  reconstructed  as  an  infinite 
series  in  the  filtered  velocity  and  its  derivatives,  with  the  lowest  order  term  in  the  form 
known  as  the  nonlinear  model 


dUj  dUj 
TlJ  dxk  dxk 


(12) 


Germano  (1986)  introduced  an  exponential  filter  which  has  an  inverse  in  a  form  of  a 
differential  operator.  The  inverse  allows  to  compute  the  SGS  stress  tensor  from  the 
definition  (3)  and  is  also  consistent  with  the  nonlinear  model  in  the  first  approximation. 
The  expression  (12)  was  first  investigated  as  a  SGS  model  by  Clark  et  al  (1979)  and 
subsequently  also,  among  others,  by  Liu  et  al.  (1991),  Borue  and  Orszag  (1998),  Leonard 
(1997),  Winckelmans  et  al  (2001),  and  Iliescu  et  al.  (2003).  In  particular,  Iliescu  et  al. 
(2003)  clarified  that  the  form  of  the  model  is  a  direct  consequence  of  approximating  the 
Fourier  transform  of  a  Gaussian  filter  by  its  Taylor  polynomial  expansion  in  filter  width 
A  and  keeping  terms  up  to  the  order  A2  in  the  expression  for  the  SGS  stress. 

The  modern  approach  to  deconvolution  methods  in  SGS  modeling,  including  the  ‘de- 
convolution’  terminology,  has  been  introduced  by  Stolz  and  Adams  (1999)  (detailed 
review  of  that  work  is  also  provided  by  Domaradzki  and  Adams  (2002)).  These  authors 
recognized  that  the  conditions  for  the  exact  inverse  of  filtering  operations  are  gener¬ 
ally  not  satisfied,  but  it  is  possible  to  use  an  approximate  or  regularized  deconvolution. 
Specifically,  Stolz  and  Adams  (1999),  and  Stolz  et  al.  (2001a)  employ  a  power  series 
expansion  of  van  Cittert  (1931),  which  assumes  that  the  inverse  filtering  operator,  if 
it  existed,  could  be  written  as  G~l  =  *“  G)n,  where  /  is  the  identity  operator. 

Regularized  inversion  is  obtained  if  this  series  is  truncated  at  a  certain  N  which  becomes 
the  regularization  parameter,  giving  an  approximate  inversion  as 


N 

Ui  %  Qn  *ut  =  Yjl  -  G)n  *  7q, 

n= 0 


(13) 


where  the  asterisk  symbolizes  the  convolution.  For  a  given  filter  the  approximate  in¬ 
verse  is  easy  to  implement  because  it  is  obtained  by  a  multiple  application  of  the  filter 
operation  to  the  filtered  function.  When  used  in  SGS  modeling  the  inverse  (13)  leads 
to  Approximate  Deconvolution  Model  (ADM)  of  Stolz  and  Adams  (1999).  A  similar 
expansion  was  used  by  Horiuti  (2001),  leading  to  a  multi-level  filtered  model. 

Other  approaches  to  the  problem  of  deconvolution  in  SGS  modeling  have  been  utilized 
though  often  without  introducing  explicitly  an  inversion  operation  (what  may  be  called 
an  ‘implicit  deconvolution’).  Shall  and  Ferziger  (1995)  expand  the  unfiltered  quantity 
in  the  integral  (4)  in  the  Taylor  series  around  an  LES  mesh  point,  which  leads  to  an 
expansion  for  the  filtered  quantity  in  terms  of  the  derivatives  of  the  exact  quantity. 
These  derivatives  can  be  subsequently  approximated  using  finite  difference  formulas  on 
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the  chosen  discretization  stencil,  leading  to  a  tridiagonal  system  of  equations  for  the 
unfiltered  quantity.  The  solution  of  the  system  provides  the  deconvolved  velocity  which  is 
then  used  in  the  definition  of  the  SGS  stress  (3).  Geurts  (1997)  constructs  an  approximate 
inverse  operator  by  requiring  that  polynomials  up  to  a  certain  order  are  recovered  exactly 
from  their  filtered  counterparts.  Domaradzki  and  Saiki  (1997)  and  Kuerten  et  al  (1999) 
work  with  a  Fourier  representation  of  a  grid  function,  i.e.,  with  discretized  version  of  (10). 
For  a  number  of  simple  filters,  such  as  the  top-hat  or  the  Gaussian,  the  inverse  for  a  single 
Fourier  mode  can  be  easily  obtained,  and  from  this  the  numerical  inverse  for  the  entire 
filtered  quantity. 

The  deconvolution  offers  a  deceptively  simple  solution  to  the  SGS  modeling  problem 
because  it  relies  only  on  the  definitions  of  the  filter  (4)  and  the  SGS  stress  tensor  (3) 
without  any  reference  to  the  physics  of  turbulence.  Since  the  need  for  the  physical  models 
is  obviated  this  approach  appears  to  be  more  promising  than  the  classical  models.  The 
key  word,  however,  is  ‘deceptively'  because  the  closer  inspection  of  the  deconvolution 
concept  reveals  that  its  simplicity  is  only  illusory.  This  is  because  the  formal  inversion 
formalism  implicitly  assumes  infinite  spectral  support  for  both  the  unfiltered  and  filtered 
quantities.  In  such  a  case,  if  an  inverse  filter  exists,  the  velocity  U{  can  be  recovered  exactly 
from  the  filtered  velocity  xit  and  Eq.  (1)  can  be  viewed  as  a  result  of  a  simple  change  of 
dependent  variables  in  the  full  Navier-Stokes  equation  from  U{  to  TL  and  both  equations 
are  formally  equivalent.  This  fact  has  long  been  recognized,  e.g.  by  Zhou  et  al.  (1989), 
Domaradzki  and  Loh  (1999),  Langford  and  Moser  (2001),  Winckelmans  et  al.  (2001),  and 
Domaradzki  and  Adams  (2002),  who  provided  critiques  of  such  exact  inversion  procedures 
as  a  SGS  modeling  tool.  The  essence  of  the  criticism  is  that  the  formal  approach  does 
not  account  for  the  fact  that  in  an  actual  LES  Eq.  (1)  is  solved  in  a  discretized  form  with 
the  numerical  resolution  much  less  than  required  to  represent  full  Navier-Stokes  solution 
ut.  This  is  because,  by  design,  the  filter  attenuates  scales  smaller  than  the  filter  width  A 
and  consequently  the  filtered  velocity  TZt  can  be  accurately  represented  on  a  mesh  with 
the  mesh  size  Ales  ~  0(  A),  which  can  be  much  larger  than  the  mesh  size  A  pus  needed 
to  represent  the  full  Navier-Stokes  field  U{  in  a  direct  numerical  simulation.  Therefore, 
in  practical  LES  in  addition  to  the  explicit  filter  (4)  an  implicit  projection  of  the  field 
Hi  is  present  that  makes  it  possible  to  represent  it  on  a  finite  mesh  with  the  mesh  size 
Ales  >  Ap^s-  In  Fourier  methods  the  LES  projection  is  equivalent  to  an  implicit  sharp 
spectral  filter  with  the  cutoff  wavenumber  kc  =  it/ Ales-  Because  of  the  implicit  filter 
the  spectral  support  for  the  LES  velocity  ut  is  much  less  than  for  the  full  velocity  Ui  and 
the  formal  inversion  of  the  filtering  operation  can  only  recover  the  full  velocity  only  up 
to  the  same  truncation  wavenumber  kc  as  the  original  LES  field  ul. 

The  traditional  notation  in  LES  usually  combines  the  effects  of  the  explicit  spatial 
filtering  and  of  the  implicit  projection  onto  a  grid  into  one  symbol,  the  overbar,  creating 
a  potential  for  misunderstanding  of  the  deconvolution  formalism  in  the  SGS  modeling. 
In  discussing  the  deconvolution  approaches  it  is  thus  beneficial  to  separate  both  opera¬ 
tions  using  an  explicit  notation  as  proposed,  for  instance,  by  Winckelmans  et  al.  (2001), 
Carat i  et  al.  (2001),  and  Domaradzki  and  Adams  (2002).  The  latter  reference  reserves 
the  overbar  for  the  explicit,  spatial  filtering  and  the  spectral  LES  truncation  is  denoted 
by  the  superscript  and  its  complement  by  ‘ S\  i.e.,  the  full  velocity  field  is  decomposed 
as 


Ui  —  uf  +  uf  . 


(14) 


Assuming  that  the  spatial  filtering  and  the  spectral  truncation  commute,  both  opera- 
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tions  are  applied  to  the  Navier-Stokes  equation  giving 


_ r  yj  7 - r  £  1  d  r 

Ul  +  (u-uj)  =  +  U- 


d2 


dt 


Ox 


Ulc 


pdxf  dXjdXj  1  * 


(15) 


The  filtered  and  truncated  product  of  velocities  in  (15)  can  be  rewritten  as  follows 

rC 


- - r£  / _ c _ C\£ 

(uxUj)  =  [Ui  Uj  ) 

giving  the  following  form  of  the  LES  equation 


(ufuf )  -  ( UiC  UjC)C  4-  [ui’uf  4-  ufuf  4-  uf  uf  j  ,  (16) 


iu<c  +  (U‘C  uic)c  =  +  vizhzruf  -  +  rdrp)' 


d2 


d 


at. 


dxj 


pdxx 


dxjdxj  dxj 


where 


and 


r^  =  (ufUf)  -(W£HJ£)  , 


(17) 


(18) 


TiJP  ~  (Uf  ^ uf  +  +  . 


(19) 


The  above  form  of  the  LES  equations  shows  that  the  SGS  stress  tensor  consists  of  two 
parts:  which  results  from  the  nonlinear  interactions  among  scales  represented  on  a 

numerical  grid  and  r™rp ,  which  accounts  for  the  the  nonlinear  interactions  that  involve 
scales  that  are  not  represented  on  a  grid. 

If  the  spatial  filtering  operation  is  invertible,  uf  can  be  obtained  from  xTxc,  and  r^p 
can  be  evaluated  without  any  SGS  modeling.  However,  in  such  a  case  the  term  rf*rp 
vanishes  because  defiltering  does  not  recover  subgrid  scales  S  lost  due  to  the  projection 
of  all  quantities  onto  a  coarse  LES  grid.  Therefore,  r[jP  would  be  the  total  SGS  stress 
that  can  be  obtained  from  a  known  LES  velocity  using  the  exact  deconvolution.  If  r^rp 
is  neglected  in  the  LES  equation  (15),  the  entire  equation  can  be  deconvolved  giving 


y  +  ±.(ufnff 


dxx 


1  d  c  d2 
pdx ^  dxjdxj 


(20) 


which  is  the  Navier-Stokes  equation  for  the  resolved  velocity  uf .  Therefore  the  application 
of  such  an  exact  deconvolution  in  LES  is  no  different  from  solving  the  equivalent,  Navier- 
Stokes  equation  for  the  de-filtered  velocity  on  a  domain  truncated  to  the  LES  cutoff.  In 
terms  of  numerical  simulations  it  means  that  LES  is  equivalent  to  DNS  performed  on 
the  coarse  LES  mesh  and  the  LES  results  can  be  simply  obtained  from  the  coarse-mesh 
DNS  data  by  the  spatial  filtering.  At  higher  Reynolds  numbers  such  DNS  are  generally 
under-resolved.  In  the  case  of  isotropic  turbulence  at  very  high  Reynolds  numbers  the 
dynamics  of  the  resolved  modes  in  such  simulations  will  be  essentially  inviscid,  leading 
to  the  cquipartition  of  energy.  In  such  a  state  the  energy  flux  between  different  modes 
vanishes  on  average  and  so  does  the  associated  SGS  dissipation  in  the  equivalent  LES. 
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Therefore,  the  pure  deconvolution  formalism,  when  implemented  on  a  coarse  LES  mesh 
will  fail  as  a  subgrid  scale  modeling  tool  for  high  Reynolds  number  flows. 

The  main  reason  why  the  deconvolution  models  in  actual  numerical  implementations 
seem  to  work  is  because  they  do  not  evaluate  rf*7*  exactly  but  rather  approximate  it,  in 
the  process  generating  an  implicit  expression  of  the  similar  form  as  (19).  For  instance, 
the  approximate  deconvolution  model  (ADM)  of  Stolz  and  Adams  (1999)  in  its  first 
formulation  did  not  provide  an  explicit  model  for  the  term  r™p  and  yet  excellent  results 
were  obtained  for  decaying  compressible  isotropic  turbulence.  In  ADM  an  approximation 
for  is  obtained  from  (13) 


U{ 


N 


rX  _ 


=  ]P(/  ~  G)n  *  ui  -  Qn  *  uf 


n=0 


(21) 


where  tilde  in  uf  indicates  the  approximation  to  the  actual  large  scale  component  of 
the  velocity  uf  that  would  be  obtained  from  the  exact  deconvolution.  The  leading-order 
term  of  the  deconvolution  error  is 

or(iV+l),,£ 

inf  =  if  -  uf  =  (22) 

dxj  ' 

where  Cj  depends  on  the  filter  kernel  which  is  of  order  r,  and  a  tensorial  extension  of  (4) 
to  multiple  dimensions  is  assumed  (Domaradzki  and  Adams  2002).  Inserting  uf  for  uf 
into  (18)  results  in 


£ 

TijP  =  (uiuj  )C  ~  (^7£  +  ( ufSuf  +  Sufuf  +  SufSuf'j  (23) 

with  the  last  term,  which  is  due  to  the  regularized  deconvolution,  providing  a  model 
for  rt"rp.  Clearly,  such  a  model  will  depend  on  the  primary  filter  and  details  of  the 
deconvolution,  e  g.,  parameter  AT,  and  thus  lacks  a  physical  basis.  In  high  Reynolds 
number  flows  ADM  suffers  from  insufficient  dissipation  and  must  be  supplemented  by 
a  dissipative  model  for  the  term  r™rp  in  (17),  a  problem  similar  to  that  faced  by  the 
classical  similarity  models.  ADM  with  the  relaxation  term  designed  to  provide  additional 
dissipation  has  been  proposed  and  evaluated  by  Stolz  et  al.  (200 1  a, b). 

Thus  the  primary  conclusion  concerning  defiltering  procedures  in  LES  is  that  such 
LES  are  fundamentally  equivalent  to  under-resolved  DNS  because  they  neglect  the  term 
r^rp  in  Eq.  (17)  which  is  the  only  term  that  provides  a  physically  realistic  mechanism 
for  the  SGS  dissipation.  This  term  accounts  for  interactions  with  the  subgrid  scales 
with  sizes  below  those  that  can  be  represented  on  the  LES  mesh  and  which  cannot 
be  recovered  by  the  deconvolution  procedures.  In  practical  LES  this  equivalence  may 
be  obscured  by  additional  numerical  effects,  e.g.,  an  approximate  deconvolution  may 
provide  a  dissipation  mechanism  absent  for  the  exact  deconvolution.  Nevertheless  the 
subgrid  scales  not  represented  on  the  numerical  grid  cannot  be  produced  by  the  formal 
inversion  of  a  filter,  either  exact  or  approximate  and,  in  general,  there  is  no  substitute 
for  modeling  their  effects.  In  that  respect  the  deconvolution  models  are  no  different  from 
the  similarity  models:  both  do  not  contain  information  about  actual  subgrid  scales  and 
because  of  that  require  extraneous  dissipation  mechanisms,  provided  by  eddy  viscosity 
or  relaxation  expressions. 
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A  different  application  of  an  explicit,  convolution  filter  and  its  approximate  inversion 
is  regularization  of  the  Navier-Stokes  equations  obtained  by  filtering  not  the  full  equa¬ 
tions  but  only  the  advective  velocity,  known  as  the  alpha-model  (Foias  et  al.  (2001)). 
Holm  and  Domaradzki  (2001)  and  Geurts  and  Holm  (2003)  showed  that  the  regularized 
Navier-Stokes  equations  can  be  transformed  into  the  usual  LES  form  with  the  specific 
expressions  for  the  SGS  stress  tensor.  These  expression  are  reminiscent  of  the  generalized 
similarity /deconvolution  models  and  thus  both  types  of  models  are  expected  to  behave 
in  the  same  manner  in  actual  LES. 


4.  Subgrid  Scale  Estimation  Models 

In  order  to  obtain  nonvanishing  T™rp  the  estimation  procedures  provide  expressions  for 
the  the  subgrid  velocity  u ^  in  terms  of  the  large  scale  velocity  u^.  Such  a  modeled  velocity 
is  sometimes  called  the  estimated  subgrid  scale  velocity  or  the  synthetic  velocity.  The 
total  modeled  velocity  is  then  simply  a  sum  of  the  known  uj  and  the  estimated  SGS 
velocity  u $  and  is  used  to  compute  the  SGS  stress  tensor  from  its  definition  (3). 

Scotti  and  Meneveau  (1999)  construct  the  subgrid  scales  from  the  known  resolved  scales 
using  the  fractal  interpolation  technique.  The  simplest  application  of  the  procedure  gives 
the  model  that  does  not  produce  enough  SGS  dissipation.  This  deficiency  can  be  corrected 
by  tuning  model  parameters  to  get  the  best  alignment  of  the  computed  SGS  stress  with 
the  resolved  rate-of-strain.  Enforcing  the  best  alignment  is  similar  to  the  assumption 
underlying  the  eddy  viscosity  models.  Another  example  of  using  synthetic  fields  is  the 
so-called  additive  turbulent  decomposition  proposed  by  Hylin  and  McDonough  (1999). 
Separate  equations  are  written  for  large  and  small  scales  in  the  decomposition  and  the 
equations  for  the  large  scales  depend  explicitly  on  the  small  scales.  The  small,  synthetic 
scales  are  constructed  as  stochastic  variables  determined  from  one-dimensional  chaotic 
maps,  substituted  into  equations  for  the  large  scales,  and  the  equations  are  advanced  in 
time.  Kerstein  (1999)  developed  the  one-dimensional  turbulence  (ODT)  model  which  has 
been  evaluated  in  LES  of  free  shear  flows  and  wall  bounded  flows  (Kerstein  et  al.  2001, 
Kemenov  and  Menon  2002).  The  essence  of  ODT  is  to  represent  a  three-dimensional 
turbulent  field  through  a  one  dimensional  line  of  data,  which  in  LES  extends  across  a 
volume  determined  by  the  neighboring  LES  niesli  points.  The  flow  variables  are  defined 
along  the  line  with  the  numerical  resolution  sufficient  to  explicitly  account  for  viscous 
effects  and  provide  a  model  for  SGS  quantities.  The  modeled  SGS  quantities  are  then 
used  to  calculate  SGS  fluxes  needed  to  advance  in  time  the  LES  equations. 

Rather  than  using  one  dimensional  procedures,  other  approaches  work  with  three  di¬ 
mensional  models  of  subgrid  scales.  Pullin  and  Saffman  (1994)  proposed  that  subgrid 
scales  are  a  collection  of  vortical  structures  with  prescribed  properties.  Stretching  of  these 
vortices  by  the  resolved  field  provides  the  mechanism  for  the  energy  transfer  from  the 
resolved  to  subgrid  scales,  i.e.,  the  mechanism  for  the  SGS  dissipation.  With  additional 
assumptions  the  analytical  formulas  for  the  SGS  vortices  allow  to  derive  an  expression 
for  the  SGS  stress  tensor  which  can  be  used  in  the  LES  equations.  The  model  has  been 
implemented  in  the  spectral  space  by  Misra  and  Pullin  (1997)  and  in  the  physical  space 
representation  by  Voelkl  et  al.  (2000),  and  also  applied  to  a  passive  scalar  mixing  by 
Pullin  (2000).  Modeling  of  subgrid-scale  vorticity  has  been  also  explored  by  Kerr  and 
Domaradzki  (1996).  In  that  work  the  subgrid-scale  vorticity  is  confined  to  a  range  of 
wavenumbers  outside  the  resolved  range,  used  for  dc-aliasing  of  pseudo-spectral  simu¬ 
lations.  The  model  is  based  on  the  vorticity  production  in  that  range  by  the  resolved 


October  20,  2010  15U  International  Journal  of  Computational  Fluid  Dynamics  Do¬ 

rn  ar  ad  zki°  re  view' ijcfd'R2 

Large  Eddy  Simulations  without  Explicit  Eddy  Viscosity  Models  11 


scales.  The  modeled  subgrid-scale  vorticity  was  used  to  compute  the  SGS  quantities 
which  compared  favorably  with  the  exact  quantities  in  a  pnon  tests. 

The  subgrid- scale  estimation  model  introduced  by  Domaradzki  and  Saiki  (1997)  at¬ 
tempts  to  model,  or  estimate,  the  unfiltered  velocity  field  appearing  in  the  definition  of 
the  subgrid-scale  stress  tensor.  It  consists  of  two  steps.  In  the  kinematic  step  an  approx¬ 
imate  inversion  of  the  filtering  operation  is  performed  producing  defiltered  field  uf.  The 
next,  dynamic,  or  nonlinear,  step,  generates  the  subgrid  scales  twice  smaller  than  the 
smallest  scales  resolved  on  the  LES  mesh,  producing  a  perturbation  velocity  uf .  In  order 
to  be  able  to  represent  these  subgrid  scales  a  mesh  finer  than  the  LES  mesh  is  explicitly 
introduced  by  halving  the  coarse  LES  mesh  in  each  Cartesian  direction.  The  full  velocity 
U{  is  approximated  by  the  estimated  velocity  ux 


ux  %  ux  —  uf  +  uf.  (24) 

The  estimated  velocity  on  the  fine  mesh  is  used  to  compute  the  SGS  stress  directly 
from  the  definition  on  the  coarser  LES  mesh.  The  SGS  estimation  model  was  applied  to 
LES  of  low  Reynolds  number  incompressible  turbulence  in  channel  flow  by  Domaradzki 
and  Saiki  (1997),  Domaradzki  and  Loh  (1999),  Loh  and  Domaradzki  (1999),  to  spatially 
evolving  compressible  turbulence  by  Domaradzki  et  al.  (1998),  Dubois  et  al  (2002),  and 
to  Rayleigh-Benard  convection  by  Kimmel  and  Domaradzki  (2000). 

The  models  described  above  share  frequently  with  the  similarity  and  deconvolution 
models  the  same  problem  of  the  insufficient  SGS  dissipation.  The  reason  for  that,  can  be 
understood  by  noting  that  the  SGS  dissipation  reflects  coupling  between  resolved  and 
subgrid  scales  which  is  induced  by  the  Navier-Stokes  dynamics.  For  instance,  it  is  known 
that  velocity  fields  obtained  in  spectral  DNS  lose  the  ability  to  sustain  the  energy  flux 
characteristic  of  turbulence  if  phases  of  the  complex  velocity  modes  are  randomized,  even 
if  their  amplitudes  are  preserved.  Only  the  dynamics  through  Navier-Stokes  equations, 
acting  over  sufficiently  long  time,  will  induce  the  proper  phase  relationships  between  re¬ 
solved  and  subgrid  scales  that  are  required  for  the  turbulent  energy  cascade.  Therefore, 
the  insufficient  SGS  dissipation  in  the  estimation  models  implies  that  the  estimated  sub¬ 
grid  scale  quantities  are  deficient  in  modeling  the  phase  relationships  that  are  consistent 
with  the  Navier-Stokes  dynamics.  In  the  next  section  we  describe  models  that  address 
this  deficiency. 


5.  Explicit  Filtering,  Implicit  Dissipation,  and  Constrained  SGS 
Models 

Domaradzki  and  Yee  (2000)  proposed  an  approach  that  addresses  the  insufficient  SGS 
dissipation  in  the  original  estimation  model  at  high  Reynolds  numbers  without  intro¬ 
ducing  eddy  viscosity  expressions.  The  method  is  based  on  the  observation  that  in  the 
original  estimation  procedure  a  single  act  of  nonlinear  interactions  among  large  scale, 
obtained  through  a  single  time  step  computation  of  the  nonlinear  term,  determines  the 
modeled  perturbation  velocity  uf.  However,  the  actual  process  of  generation  of  sub¬ 
grid  scales  involves  much  larger  number  of  nonlinear  interactions,  i.e.  involving  multiple 
computations  of  nonlinear  term  over  many  time  steps,  covering  approximately  one  large 
eddy  turnover  time.  This  suggests  that  the  realism  of  the  modeled  velocity  field  can  be 
improved  by  allowing  additional  nonlinear  interactions  in  its  generation.  This  goal  is 
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accomplished  by  allowing  the  total  estimated  velocity  uu  Eq.  (24)  to  evolve  according  to 
the  truncated  Xavier-Stokes  (TNS)  equations  on  the  fine  mesh 


5*  + 


1  a  _  a2  _ 

- V  i-  V - U‘ 

pdii  dxjdxj 


(25) 


—u,  =  0.  (26) 

The  term  ‘truncated’  applied  to  equations  (25)-(26)  refers  not  to  the  form  of  the  equa¬ 
tions,  which  is  no  different  from  the  form  of  the  full  Xavier-Stokes  equations,  but  to 
the  limited  numerical  resolution  used  to  solve  these  equations.  The  procedure  consists  of 
several  steps,  symbolically  represented  as  follows 

Ui(t)  der-%lv-  [uf  (t)  +  uf  (<)]  N-*S  u, (t  +T)  -¥  [ u,(t  +  T)  -  uf  (<)]  Ui(t  +  T). 

(27) 

The  LES  velocity  U{  at  time  t  is  given  on  a  coarse  mesh  and  is  first  defiltered  through  an 
exact  or  approximate  deconvolution,  giving  uf .  This  field  is  interpolated  to  a  finer  mesh, 
obtained  normally  by  halving  a  mesh  size  for  the  coarse  mesh,  and  the  SGS  velocity  uf  is 
generated  using  the  velocity  estimation  model.  The  combined  field  is  then  advanced  over 
time  T  using  truncated  Xavier-Stokes  equations  (25)-(26),  which  results  in  enhancing 
nonlinear  couplings  between  large  and  small  scales.  The  purpose  of  the  remaining  two 
steps  is  to  obtain  the  filtered  LES  velocity  after  the  time  period  T.  Note  that  rather  than 
simply  filtering  the  Xavier-Stokes  solution  ut(t  T),  the  original  perturbation  velocity 
at  time  t,  uf{t),  is  first  subtracted,  followed  by  filtering  and  sampling  of  the  result 
on  the  coarse  mesh.  The  subtraction  is  introduced  to  avoid  spurious  dynamics  by  the 
combination  of  the  estimation  step  and  the  reduction  step.  This  way  if  the  procedure  is 
applied  without  the  Xavier-Stokes  evolution  the  initial  field  ut  is  recovered  exactly  at  the 
end  of  the  reduction  step,  i.e.,  modifications  of  the  filtered  field  TZj  axe  caused  only  by  the 
truncated  Xavier-Stokes  dynamics  of  the  estimated  velocity  Ui.  For  flows  at  high  Reynolds 
numbers  the  fine  mesh  used  is  still  far  too  coarse  to  resolve  all  relevant  scales  and  the 
truncated  Xavier- Stokes  equations  would  eventually  lead  to  the  unphysical  equipartition 
of  energy  at  large  wavenumbers.  The  filtering  and  sampling  to  a  coarse  LES  mesh,  applied 
every  period  T,  prevents  the  accumulation  of  the  energy  in  the  small  scales.  Therefore, 
the  time  T  must  be  short  enough  to  prevent  equipartition  but  long  enough  to  allow  the 
sufficient  nonlinear  couplings  between  large  and  small  scales  to  develop.  It  was  found 
that  the  time  T  can  be  chosen  as  the  eddy  turnover  time  for  the  smallest  resolved  scales. 
Because  the  Xavier-Stokes  equations  are  used  to  induce  nonlinear  couplings  between 
the  large  and  small  scales  and  thus  the  corresponding  energy  flux,  the  method  is  less 
dependent  on  the  accuracy  of  the  prediction  for  the  initial  estimated  velocity  uf(t)- 
Even  setting  uf  (t)  simply  to  zero  and  allowing  the  truncated  Xavier-Stokes  equations  to 
generate  subgrid  scales  over  time  T  can  often  be  sufficient  to  obtain  accurate  LES  results 
as  shown  by  Domaradzki  and  Horiuti  (2001)  and  Yang  and  Domaradzki  (2004). 

Operationally,  the  LES  procedure  (27)  consists  of  a  number  of  Xavier-Stokes  runs 
performed  on  the  fine  mesh  without  any  explicit  model  and  re-initializatiou,  involving 
filtering  and  sampling  to  the  coarser  LES  mesh,  every  time  period  T.  Since  all  quantities 
between  defiltering  and  filtering  steps  are  represented  on  the  fine  mesh,  and  the  Xavier- 
Stokes  equations  are  solved  on  the  same  mesh,  the  method  can  be  further  simplified  by 
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removing  the  deconvolution  part  of  the  procedure  entirely  as  suggested  by  Domaradzki 
et  al.  (2002).  This  modification  uses  only  the  large  scale  velocity  obtained  from  the 
TNS  solution  by  the  appropriate  filtering,  symbolically 


►[uf  {t)  4*  u]p(£)]  Ui(t  +  T)  ^ — >  u^it  +  T). 


(28) 


In  simulations  using  (28)  only  a  single,  fine  mesh  is  employed.  For  high  Reynolds  number 
flows  the  TNS  method  is  thus  equivalent  to  a  sequence  of  under-resolved  Navier  Stokes 
simulations  with  periodic  filtering  of  the  numerical  solution.  There  are  several  advan¬ 
tages  of  such  an  approach:  since  only  a  solution  is  periodically  modified,  an  underlying 
Navier-Stokes  solver  can  be  used  in  the  same  form  for  DNS  and  LES;  the  method  is 
easily  generalized  to  model  other  nonlinear  phenomena,  e.g.,  convective  or  compressible 
turbulent  flows,  by  periodically  filtering  corresponding  flow  variables;  the  difficulties  with 
transformation  properties  for  the  modeled  SGS  stresses  (Horiuti  2001)  are  avoided.  The 
obvious  problem  with  this  approach  is  that  the  periodic  filtering  generates  a  discontinu¬ 
ous  in  time  solution  Uj.  This  problem  is  less  serious  than  it  appears  because  in  LES  one 
is  interested  only  in  the  large  scale  component  ,  while  the  small  scale  component  uf 
is  considered  to  be  a  model;  by  design,  typical  filters  used  affect  only  the  small  scales  in 
the  method  and  discontinuity  in  the  large  scales  is  minimized.  The  discontinuity  can  be 
entirely  removed  if  the  procedure  is  formulated  in  the  framework  of  the  LES  equations 
coupled  with  the  TNS  (Domaradzki  et  al.  2002)).  The  TNS  equations  are  advanced  in 
time  and  at  each  time  step  a  SGS  stress  tensor  is  calculated  using  the  velocity  ux 
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where  the  notation  indicates  that  the  SGS  stress  is  computed  from  the  estimated  velocity 
iq  given  on  the  fine  mesh,  using  the  overbar  filter  and  sampling  on  a  coarser  LES  mesh. 
The  SGS  stress  computed  this  way  is  used  as  the  approximation  to  the  actual  SGS  stress 
to  advance  in  time  the  LES  equations  (1)  discretized  on  the  LES  inesh.  The  discontinuity 
in  Ui,  introduced  by  the  periodic  in  time  filtering,  is  considerably  reduced  by  the  overbar 
filtering  and  sampling  involved  in  computing  (29).  This  results  in  a  smooth  evolution  of 
the  LES  velocity  Uj,  more  so  that  the  SGS  term  in  LES  equations  is  much  smaller  than 
the  resolved  nonlinear  term.  Such  considerations  imply  that  the  implementation  of  TNS 
described  above  is  preferred;  however,  it  is  quite  cumbersome  and  costly  because  two 
sets  of  equations  on  two  different  meshes  must  be  solved  numerically.  Because  of  that, 
in  practice  the  implementation  (28)  on  a  single,  fine  mesh  is  used  and  the  discontinuity 
in  the  numerical  solution  is  tolerated. 

The  TNS  method  depends  formally  on  the  filter  used  and  the  filtering  interval  T.  The 
filter  dependence  has  not  been  investigated  in  detail  but  the  dependence  appears  weak 
for  filters  that  affect  only  the  small  scales,  determined  by  'wavenumbers  between  1/2 kc 
and  fcc,  where  kc  is  the  mesh  cutoff.  The  filtering  interval  is  a  free  parameter,  normally 
chosen  by  trial  and  error,  using  the  physical  guideline  that  it  is  determined  by  the  eddy 
turnover  time  of  the  subgrid  scales.  Once  the  optimal  selection  is  made,  the  results  of 
LES  are  only  weakly  dependent  on  the  filtering  interval;  the  variation  within  a  factor 
of  two  gives  acceptable  results.  Recently,  Tantikul  and  Domaradzki  (2010)  have  shown 
that  for  a  prescribed  filter,  the  filtering  interval  T  can  be  determined  automatically  in 
the  course  of  the  simulations  using  a  criterion  based  on  the  energetics  of  the  modeled 
subgrid  scales.  This  way,  for  a  given  filter,  the  method  becomes  parameter-free. 
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Another  way  of  modifying  the  estimation  concept  to  produce  sufficient  SGS  dissipa¬ 
tion  was  used  by  Park  and  Mahesh  (2008).  They  proposed  the  resolved  subgrid  scale 
estimation  model  (RSEM)  which  constraints  the  SGS  stress  to  match,  in  the  mean  or 
least  squares  sense,  the  SGS  dissipation  of  the  dynamic  Smagorinsky  model  for  the  same 
velocity  field.  The  model  does  not- have  the  eddy  viscosity  form,  but  the  extraneous  con¬ 
straint  guarantees  that  the  total  SGS  dissipation  is  comparable  to  that  produced  by  a 
selected  eddy  viscosity  model. 

The  fact  that  TNS  solves  Navier-Stokes  equations  without  any  modeling  terms  is  su¬ 
perficially  similar  to  the  so-called  implicit  LES,  the  approach  that  was  originally  proposed 
by  Boris  et  al.  (1992)  and  reviewed  recently  in  the  monograph  edited  by  Grinstein  et  al. 
(2007).  In  ILES  the  Navier-Stokes  equations,  also  without  any  SGS  models,  are  solved 
numerically  on  a  coarse  LES  mesh  using  higher-order  nonoscillatory  methods,  such  as 
total  variation  diminishing  (TVD),  flux-correct ed-transport  (FCT)  and  flux-limited  and 
sign-preserving  schemes  (Sweby  1984,  Zalesak  1979,  Harten  et  al.  1987)),  originally  de¬ 
veloped  to  control  numerical  oscillations  in  problems  involving  steep  gradients  and/or 
discontinuities  such  as  shocks.  The  ILES  methodology  was  initially  justified  on  the  basis 
of  a  practical  observation  that  truncation  errors  in  such  discretizations  of  Navier-Stokes 
equations  introduce  numerical  dissipation  and  its  effect  is  qualitatively  similar  to  the 
effects  of  the  explicit  SGS  models.  For  instance,  Porter  and  Woodward  (1994)  report  the 
development  of  the  fc“5/3  energy  spectrum  in  numerical  simulations  of  decaying  isotropic 
turbulence  performed  using  piecewise  parabolic  method  (PPM)  implemented  in  an  Euler 
solver,  i.e.,  for  equations  without  any  explicit  viscous  terms.  While  such  results  can  be 
used  to  support  the  ILES  methodology,  it  is  also  possible  to  find  counterexamples.  For 
instance,  Gamier  et  al.  (1999)  analyzed  several  different  shock-capturing  Euler  schemes 
applied  to  decaying  isotropic  turbulence  and  the  conclusions  were  less  favorable.  While 
it  was  possible  to  obtain  the  inertial  subrange,  the  other  results,  e.g,  pdfs  of  the  velocity 
derivatives  and  pressure  showed  a  behavior  typical  of  low  Reynolds  number  flows  rather 
than  expected  from  high  Reynolds  number  LES.  This  behavior  was  attributed  to  the 
fact  that  the  numerical  dissipation  was  often  much  greater  than  the  SGS  dissipation 
computed  for  the  same  field  using  the  Smagorinsky  model.  Based  on  these  results  ILES 
may  provide  substantially  more  numerical  dissipation  than  expected  based  on  the  physics 
of  turbulent  cascade.  Similar  conclusions  were  reached  by  Domaradzki  and  Radhakrish- 
nan  (2005)  who  showed  that  the  ILES  results  for  rotating  and  non-rotating  turbulence 
were  sensitive  to  the  time  step  and  the  method  failed  to  produce  theoretically  expected 
results  for  certain  initial  conditions  and  for  rotating  turbulence.  Therefore,  while  ILES 
offers  the  ease  of  implementation  because  only  Navier-Stokes  equations  are  being  solved, 
the  uncertainty  regarding  the  numerical  dissipation  may  negate  potential  advantages. 

Recently,  Hickel  et  al.  (2006)  demonstrated  that  the  ILES  methodology  can  be  made 
viable  by  developing  a  specific  method  that  respects  the  physical  requirements  for  the 
energy  transfer  in  turbulence.  The  method  is  based  on  a  new  nonlinear  discretization 
scheme,  the  adaptive  local  deconvolution  method  (ALDM),  which  contains  several  free 
deconvolution  parameters  that  allow  to  control  the  truncation  error  of  the  scheme.  The 
free  parameters  are  constrained  such  that  the  numerical  viscosity  optimally  matches  the 
theoretical  eddy  viscosity  predicted  by  the  analytical  theories  of  turbulence.  While  the 
optimization  was  performed  for  isotropic  turbulence,  the  parameters  of  the  scheme  can 
be  used  to  simulate  other  turbulent  flows  without  explicit  SGS  models. 

ILES  and  TNS  are  related  to  other  procedures  for  underresolved  simulations  of  tur¬ 
bulence  that  rely  on  purely  numerical  techniques  to  achieve  stable  simulations.  In  the 
so-called  stabilized  spectral  LES  (Minguez  et  al.  2009)  the  numerical  stability  is  not  pro- 
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vided  by  the  truncation  error  of  the  numerical  discretization  (which  is  exponentially  small 
for  a  spectral  method  (Boyd  2001))  but  by  the  spectral  filter  applied  at  each  time  step 
that  strongly  attenuates  the  small  resolved  scales.  In  the  same  spirit  Bogey  and  Bailly 
(2006)  use  an  explicit  filter  applied  every  few  time  steps  as  a  substitute  for  a  SGS  model 
in  LES  of  a  turbulent  jet  flow.  However,  one  must  be  aware  that  simply  guaranteeing 
numerical  stability  does  not  guarantee  physically  correct  dynamics  of  the  resolved  scales. 
Therefore,  the  results  from  such  LES  have  to  be  always  compared  with  experiments  and 
simulations,  either  fully  resolved  DNS,  or  LES  performed  with  other  models,  to  gain 
confidence  that  the  method  is  not  only  numerically  stable  but  also  physically  correct. 

In  this  section  wc  have  focused  on  procedures  that  use  either  explicit  filters  or  im¬ 
plicit  dissipation  and/or  constraints  that  result  in  the  SGS  dissipation  of  a  model  to  be 
consistent  with  the  turbulence  dynamics.  It  appears  that  such  approaches  offer  great¬ 
est  promise  in  developing  viable  methods  that  are  sufficiently  dissipative  and  do  not 
incorporate  explicit  eddy  viscosity  expressions  in  the  governing  LES  equations. 


6.  Conclusions 

The  traditional  models  of  subgrid  scale  terms  such  as  a  stress  tensor  or  a  heat  flux  usually 
employ  eddy  viscosity  concepts.  The  primary  contribution  of  the  eddy  viscosity  models 
is  to  provide  the  mean  SGS  dissipation,  which  can  represent  accurately  the  global  energy 
flux  if  the  model  constants  are  properly  chosen.  However,  many  other  features  of  SGS 
terms  are  poorly  represented  by  the  eddy  viscosity  models,  providing  a  motivation  for 
research  on  alternative  approaches.  In  this  paper  we  have  focused  on  those  modeling 
approaches  that  specifically  exclude  any  use  of  the  explicit  eddy  viscosity  expressions  in 
the  LES  equations,  either  in  the  main  or  in  the  supporting  role.  Such  an  exclusion  allows 
to  evaluate  merits  of  each  approach  in  its  own  right.  We  have  discussed  the  similarity 
models,  the  deconvolution  and  nonlinear  models,  the  estimation  models,  the  implicit  (nu¬ 
merical)  models,  and  the  explicit  filtering  of  a  numerical  solution  There  arc  a  number  of 
advantages  offered  by  such  methods.  For  instance,  most  of  these  models  naturally  provide 
both  the  dissipative  forward  transfer  and  the  anti-dissipative  baekscatter  of  energy,  con¬ 
trary  to  the  eddy  viscosity  models  that  are  purely  dissipative.  Such  modeling  procedures 
arc  good  candidates  for  modeling  strongly  anisotropic  flows  because  often  their  deriva¬ 
tions  do  not  require  assumptions  of  local  isotropy  and  of  the  inertial  range.  Moreover, 
the  modeling  principles  can  be  extended  without  difficulty  to  model  additional  physi¬ 
cal  effects,  e.g.  the  density  variations  in  convective,  stably  stratified,  and  compressible 
turbulence.  Finally,  such  methods  also  avoid  difficulties  associated  with  the  transforma¬ 
tion  properties  between  different  frames  of  reference  encountered  by  the  classical  eddy 
viscosity  models. 

Despite  the  potential  advantages  listed  above,  many  non-eddy  viscosity  approaches  suf¬ 
fer  from  the  insufficient  subgrid  scale  dissipation,  often  leading  to  their  failure  in  actual 
LES.  When  faced  with  the  insufficient  SGS  dissipation  the  temptation  is  to  add  an  eddy 
viscosity  model.  This,  of  course,  largely  defeats  the  idea  of  developing  new  SGS  models 
because  it  is  clear  that  it  is  the  presence  of  the  eddy  viscosity  expressions  that  makes  the 
combined  models  work  in  actual  LES.  A  better  approach  might  be  to  try  to  understand 
the  fundamental  reasons  for  observed  failures  in  the  hope  that  they  can  be  corrected. 
However,  the  reasons  for  the  observed  failures  are  not  always  clear  in  the  literature  on  the 
subject.  One  of  the  goals  of  this  paper  was  to  provide  heuristic  explanations.  The  similar¬ 
ity  model  was  shown  to  be  deficient  in  removing  the  energy  from  the  vicinity  of  the  LES 
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cutoff  and  compounding  the  problem  by  increasing  the  energy  loss  by  the  large  scales.  In 
a  particular  case  of  an  exact  deconvolution  the  insufficient  SGS  dissiption  was  shown  to 
be  a  result  of  the  fundamental  equivalence  between  the  model  and  under  resolved  direct 
numerical  simulations.  The  estimation  model  includes  scales  beyond  the  resolution  of 
an  LES  mesh  which,  in  principle,  should  provide  a  physically  viable  mechanism  for  the 
SGS  energy  transfer.  However,  even  in  this  case  the  estimation  models  were  found  not  to 
provide  enough  SGS  dissipation  for  flows  at  high  Reynolds  numbers.  In  the  algorithmic 
estimation  procedures  the  source  of  insufficient  SGS  dissipation  usually  is  attributed  to 
an  inability  of  the  models  to  predict  the  proper  phase  relationship  between  known,  re¬ 
solved  scales  and  unknown  subgrid  scales.  This  relationship  is  a  result  of  a  complicated 
process,  involving  multiple  nonlinear  interactions  between  resolved  and  subgrid  scales, 
acting  over  at  least  one  large  eddy  turnover  time,  as  well  as  of  the  pressure  effects.  These 
effects  can  only  be  captured  using  the  dynamic  equations  for  the  subgrid  scales  as  is 
done  in  the  TNS  approach.  It  is  also  possible  to  design  an  estimation-like  model  with  the 
sufficiently  high  SGS  dissipation  if  an  extraneous  condition  is  used  that  constrains  the 
SGS  dissipation  of  the  model  by  the  SGS  dissipation  of  the  eddy  viscosity  model  for  the 
same  velocity  field.  ILES  is  another  class  of  models  that  do  not  use  explicit  eddy  viscosity 
expressions.  In  ILES  only  Navier-Stokes  equations  are  being  solved,  thus  offering  a  rela¬ 
tive  ease  of  the  implementation,  but  this  advantage  is  largely  negated  by  the  uncertainty 
how  the  numerical  dissipation  is  related  to  the  actual  SGS  energy  flux.  The  promising 
way  around  this  problem  is  to  constrain  the  numerical  dissipation  to  be  consistent  with 
the  eddy  viscosity  predicted  by  the  analytical  theories  of  turbulence.  Finally,  we  have 
reviewed  approaches  that  rely  on  the  explicit  filtering  of  a  numerical  solution.  They  fall 
into  two  classes.  One  uses  explicit  filtering  to  simply  stabilize  the  numerics  and  counts  on 
the  correct  energy  removal  rate  from  the  system  as  a  byproduct,  making  such  approaches 
similar  to  ILES.  Another,  employs  filtering  coupled  with  the  physical  considerations  of 
the  energy  transfer  in  turbulence  to  determine  the  filter  strength  and/or  frequency  of  its 
application.  The  example  of  such  an  approach  is  TNS. 

In  this  review  we  have  focused  on  general  concepts  of  non-eddy  viscosity  approaches  to 
subgrid-scale  modeling  without  discussing  results  of  specific  LES  performed  with  these 
SGS  models.  Such  results  are  available  in  the  quoted  references  and  mostly  deal  with 
fully  turbulent  flows  though  transition  to  turbulence  is  sometimes  addressed  (Schlatter 
et  al.  (2004)).  The  non-eddy  viscosity  models  are  tested  as  any  turbulence  models 
using  comparisons  with  available  experimental  and  DNS  data.  Additionally,  since  the 
motivation  is  to  replace  eddy  viscosity  models,  comparisons  against  benchmark  LES 
results  obtained  with  eddy  viscosity  models  are  common.  The  latter  dictates  that  the 
numerical  resolution  and  the  time  step  are  the  same  or  comparable  for  both  types 
of  models.  For  instance,  tests  of  TNS  reported  recently  by  Tantikul  and  Doniaradzki 
(2010)  for  turbulent  channel  flow  for  ReT  varying  from  around  200  to  2000  use  a  grid 
that  is  appropriate  for  LES  performed  with  the  dynamic  Smagorinsky  model  (Piomelli 
(1993))  and  consistent  with  the  requirements  for  wall  resolved  LES  (Sagaut  (2002)). 
The  total  number  of  mesh  points  used  in  such  LES  can  be  as  low  as  a  few  tenths  of  one 
percent  of  a  DNS  resolution,  giving  significant  savings  in  terms  of  computational  effort. 
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